We introduce a novel numerical method for the simulation of soft granular materials, in which the particles can undergo large strains under load without rupture. The proposed approach combines an explicit total Lagrangian formulation of Material Point Method (TLMPM) with the Contact Dynamics (CD) method. The TLMPM resolves particle bulk deformations whereas the CD treats contact interactions between soft particles. The efficiency and accuracy of this approach are illustrated by analyzing diametral compression of a soft circular particle and the compaction of an assembly of soft particles up to very high levels of packing fraction. We show that although the assembly undergoes a jamming transition, the particles continue to rearrange as they get increasingly distorted under load. Interestingly, as the packing fraction increases, a transition occurs to a regime fully governed by particle shape change. The evolution of the global stress as well as the connectivity of the particles as a function of the packing fraction are discussed and a predictive model relating stress to packing fraction beyond jamming transition is proposed.