In this work, we present a new variant of the balancing domain decomposition by constraints preconditioner that is robust for multi-material problems. We start with a well-balanced subdomain partition, and based on an aggregation of elements according to their physical coefficients, we end up with a finer physics-based (PB) subdomain partition. Next, we define corners, edges, and faces for this PB partition, and select some of them to enforce subdomain continuity (primal faces/edges/corners). When the physical coefficient in each PB subdomain is constant and the set of selected primal faces/edges/corners satisfy a mild condition on the existence of acceptable paths, we can show both theoretically and numerically that the condition number does not depend on the contrast of the coefficient across subdomains. An extensive set of numerical experiments for 2D and 3D for the Poisson and linear elasticity problems is provided to support our findings. In particular, we show robustness and weak scalability of the new preconditioner variant up to 8232 cores when applied to 3D multi-material problems with the contrast of the physical coefficient up to 108108 and more than half a billion degrees of freedom. For the scalability analysis, we have exploited a highly scalable advanced inter-level overlapped implementation of the preconditioner that deals very efficiently with the coarse problem computation. The proposed preconditioner is compared against a state-of-the-art implementation of an adaptive BDDC method in PETSc for thermal and mechanical multi-material problems.