In this thesis, we discuss the development, implementation and testing of a piecewise linear (PWL) continuous Galerkin finite element method applied to the threedimensional diffusion equation. This discretization is particularly interesting because it discretizes the diffusion equation on an arbitrary polyhedral mesh. We implemented our method in the KULL software package being developed at Lawrence Livermore National Laboratory. This code previously utilized Palmer's method as its diffusion solver, which is a finite volume method that can produce an asymmetric coefficient matrix. We show that the PWL method produces a symmetric positive definite coefficient matrix that can be solved more efficiently, while retaining the accuracy and robustness of Palmer's method. Furthermore, we show that in most cases Palmer's method is actually a non-Galerkin PWL finite element method. Because the PWL method is a Galerkin finite element method, it has a firm theoretical background to draw from. We have shown that the PWL method is a well-posed discrete problem with a second-order convergence rate. We have also performed a simple mode analysis on the PWL method and Palmer's method to compare the accuracy of each method for a certain class of problems. Finally, we have run a series of numerical tests to uncover more properties of both the PWL method and Palmer's method. These numerical results indicate that the PWL method, partially due to its symmetric matrix, is able to solve large-scale diffusion problems very efficiently.