In this paper, we address the main challenges of the implementation of finite
element methods for solving spatial fractional problems on three dimensional
irregular convex regions. Different from the integer case, the non-locality of
fractional derivative operators makes the assembly of fractional stiffness
matrix much more difficult, mainly in two aspects: one is the search of the
integral path of Gaussian points, and the other is the calculation of
fractional derivatives of basis functions at Gaussian points. By introducing
the ray-simplex intersection algorithm, we present an efficient method for
finding integration paths of Gaussian points. By analyzing the expression of
the fractional derivative of finite element basis functions, we give a method
for efficiently calculating the fractional derivative. In order to speed up the
procedure in MATLAB, some implementation techniques are introduced. We
demonstrate our method by applying it to different problems, including steady
and transient spatial fractional problems.