The modeling of fracture problems within geometrically linear elasticity is often based on the space of generalized functions of bounded deformationGSBDp(Ω),p∈ (1,∞), their treatment is however hindered by the very low regularity of those functions and by the lack of appropriate density results. We construct here an approximation ofGSBDpfunctions, forp∈ (1,∞), with functions which are Lipschitz continuous away from a jump set which is a finite union of closed subsets ofC1hypersurfaces. The strains of the approximating functions converge strongly inLpto the strain of the target, and the area of their jump sets converge to the area of the target. The key idea is to use piecewise affine functions on a suitable grid, which is obtainedviathe Freudenthal partition of a cubic grid.