The major difference between convex and concave particles is that more than one contact point can be found between concave particles. The existence of a multi-point contact introduces difficulties in determining each point. Another difficulty in modeling concave particles is that the normal vectors of two concave particles at the contact point are not always parallel, thereby causing differences between the action and reaction forces. Given their concavity and local convexity, heart-shaped particles are used to simulate a concave particle system. This work proposes a grid method for calculating the coordinates of each contact point. A novel data structure of multi-point contact and a searching point algorithm are also devised. Random packings under gravity and certain normal stresses are simulated to verify the validity of this method, and periodic frictionless systems are examined to find the closest packing. The validity of this method is proven by its numerical stability and convergence, and simulations of super-ellipsoid and torus particle systems attest of the universality of this method for different-shaped particles. (C) 2020 Elsevier B.V. All rights reserved.