We present a new contact algorithm that endows the granular element method [1] with the ability to model non-convex particles using non-uniform rational basis splines. This significant extension allows for the representation of particle morphological features, namely, sphericity and angularity, to their fullest extent, with local contact rolling resistance and interlocking emanating directly from grain geometry. Both particle elasticity and friction at the contact level are treated implicitly and simultaneously, and the contact algorithm is cast into a mathematical programming-based contact dynamics framework. The framework provides the advantages of implicit time integrators (for e.g., stability and larger time steps) and ability to handle both rigid and highly stiff particles. By allowing for particle non-convexity, modeling flexibility is significantly enhanced, to a level that is comparable with isogeometric methods. As such, the transition from image data to particle shapes is greatly streamlined. More importantly, increased macroscopic strength in granular packings comprising of non-convex particles is fully captured. All the above capabilities are achieved under a very modest implementation effort.