Concrete consists of coarse aggregates, mortar matrix, and interfacial transition zone (ITZ) between them at the mesoscale. Considering these three phases, many numerical tests have been conducted to study the mesodamage behavior of concrete, in which a variety of numerical methods have also been adopted. These methods are mainly based on the finite element method (FEM); however, some other methods have been proven to be helpful as well. For example, the material point method (MPM) has the advantage of building a numerical model based on pixel or voxel of the image and is capable of solving large deformation problems. In view of this, MPM is introduced in this paper. Firstly, a method for establishing the numerical specimen is put forward, considering the original sample of its mesoscopic geometric character. Then, a stochastic damage constitutive model considering the heterogeneity of the concrete is proposed. Next, the numerical model and the constitutive model are incorporated into an MPM code to conduct numerical tests. The uniaxial tension and compression tests of a random-aggregate model and a double-aggregate specimen under uniaxial tension are then simulated numerically to validate the proposed method. Results show that the proposed method can well capture the main macroscopic mechanical behavior of concrete and the mesoscopic damage initiation and propagation. It is also found that MPM can save the time of model establishing and improve calculation efficiency. The influences of different parameters of the proposed constitutive model are also clarified through a parametric study. The proposed method can provide a useful tool for concrete numerical testing and for studying the mechanical behavior of concrete at mesoscale.