Summary We present an efficient approach to compute the matrices used in finite element formulations where self‐equilibrated (zero divergence) approximations of the stress field are used. The fundamental aspect of this approach is that it is applicable to polynomial approximations of high degree (it was tested up to degree 10), providing closed‐form expressions for the elementary matrices involved. The components of the stress field are defined as a function of the coordinates in the master element, facilitating the postprocessing of the finite element results and allowing for the consideration of geometric parameters in reduced order models.

Summary The proper generalized decomposition is a well‐established reduced order method, used to efficiently obtain approximate solutions of multi‐dimensional problems in a procedure that controls the effects of the “curse of dimensionality.” The question of assessing the quality of the solutions obtained and adapting the approximations assumed, for example, the finite element meshes used, so that the best result is obtained at minimal cost, remains a relevant challenge. This article deals with finite element solutions for solid mechanics problems, using the error obtained from a dual analysis, the difference between complementary solutions, to bound the error in the solutions and to drive an optimal adaptivity process, which obtains meshes with errors significantly lower than those obtained using a uniform refinement.