A key element in modeling and simulation of impact events is the prediction of constitutive nonlinearity for geomaterials undergoing extreme deformations that also induce widespread fracture, called rubbelization or shattering. This presentation introduces box-constrained optimization problems to describe rate-dependent bi-apex stress return algorithms (i.e., cap viscoplasticity) designed for coupling with discrete fracture models. Local constitutive solution, which relies on a modern variant of classical L-BFGS-B optimization, offers dual advantages: a well-described trade-off between accuracy and speed, and the opportunity to machine learn thermodynamically-consistent governing plastic potentials. This combination provides a straightforward workflow for the constitutive model’s prototyping and calibration but also maintains physical interpretability and computational speed of its numerical implementation. Optimization-based material modeling is then implemented in the combined finite-discrete element method (FDEM) code HOSS, in order to simulate widely-distributed energetic fracture and contact along induced surfaces coupled to ductile deformation. Numerical examples, including medium-scale fragmentation due to contact with compliant impactors, are exhibited.