In this chapter we describe the Lagrangian numerical relativity code SPHINCS_BSSN. This code evolves spacetimes in full General Relativity by integrating the BSSN equations on structured meshes with a simple dynamical mesh refinement strategy. The fluid is evolved by means of freely moving Lagrangian particles that are evolved using a modern Smooth Particle Hydrodynamics (SPH) formulation. To robustly and accurately capture shocks, our code uses artificial dissipation terms, but similar to Finite Volume schemes, we apply a slope-limited reconstruction within the dissipative terms and in addition we use time-dependent dissipation parameters, so that dissipation is only applied where needed. The technically most complicated, but absolutely crucial part of our methodology, is the coupling between the particles and the mesh. For the mapping of the energy-momentum tensor \(T_{\mu \nu }\) from the particles to the mesh, we use a sophisticated combination of “Local Regression Estimate” (LRE) method and a “multi-dimensional optimal order detection” (MOOD) approach which we describe in some detail. The mapping of the metric quantities from the grid to the particles is achieved by a quintic Hermite interpolation. Apart from giving an introduction to our numerical methods, we demonstrate the accurate working of our code by presenting a set of representative relativistic hydrodynamics tests. We begin with a relativistic shock tube test, then compare the frequencies of a fully relativistic neutron star with reference values from the literature and, finally, we present full-blown merger simulations of irrotational binary systems, one case where a central remnant survives and another where a black hole forms, and of a binary where only one of the stars is rapidly spinning.