Granular materials are made up of particles with complex non-convex shapes that in many cases become connected by cementation, sintering or other adhesive processes. This paper introduces and describes the methodology behind using the bonded particle method (BPM) with the level set discrete element method (LS-DEM-BPM), in order to model connected granular materials with arbitrary particle shape. The method is thoroughly detailed in both two and three dimensions. Examples of use for the method are shown for three distinct contexts. The first is the modeling of a fault gouge, where it is demonstrated that the method has the potential to be predictive of the rate-and-state friction law and have key insights into the micromechanics of the process. Another example is presented where a cone penetrates porous sintered ice showing that bond thickness has a considerable effect on model behavior. The third example is a simulation using many unique bonds for each contact interface, demonstrating the ability to simulate partial fractures. A discussion is also included on how and when this method can be either advantageous or disadvantageous.