Multistage manufacturing processes with identical stages provide three-dimensional process data in which the first dimension represents the process (control/sensing) variable, the second is the stage, and the third is the measurement/sampling/data acquisition time point. Diagnosing quality faults in such processes often requires the simultaneous identification of crucial process variables and stages associated with product quality anomalies. Most existing diagnosis methods convert 3D data into a 2D matrix, resulting in loss of information and reduced diagnostic accuracy and stability. To address this challenge, we propose a penalized tensor regression model that regresses the product quality index against its 3D process data. For the estimation of high-dimensional regression coefficients with the limited amount of historical data, we apply the CANDECOMP/PARAFAC and Tucker decompositions to the coefficient tensor, which significantly reduces the number of parameters to be estimated. Based on the decompositions, a new regularization term is designed to enable the joint identification of critical process variables and stages. To estimate the parameters, we develop the block coordinate proximal descent algorithm and provide its convergence guarantee. Numerical studies demonstrate that the proposed methods can enhance diagnostic stability and on average improve the diagnostic accuracy by around 20% over existing benchmarks.