Abstract A new method for numerical simulation of failure behavior, namely, FEM- β , is proposed. For a continuum model of a deformable body, FEM- β solves a boundary value problem by applying particle discretization to a displacement field; the domain is decomposed into a set of Voronoi blocks and the non-overlapping characteristic functions for the Voronoi blocks are used to discretize the displacement function. By computing average strain and average strain energy, FEM- β obtains a numerical solution of the variational problem that is transformed from the boundary value problem. In a rigorous form, FEM- β is formulated for a variational problem of displacement and stress with different particle discretization, i.e., the non-overlapping characteristic function of the Voronoi blocks and the conjugate Delaunay tessellations, respectively, are used to discretize the displacement and stress functions. While a displacement field is discretized with non-smooth functions, it is shown that a solution of FEM- β has the same accuracy as that of ordinary FEM with triangular elements. The key point of FEM- β is the ease of expressing failure as separation of two adjacent Voronoi blocks owing to the particle discretization that uses non-overlapping characteristic functions. This paper explains these features of FEM- β with results of numerical simulation of several example problems.