  In [Bonito et al., J. Comput. Phys. (2022)], a local discontinous Galerkin
method was proposed for approximating the large bending of prestrained plates,
and in [Bonito et al., IMA J. Numer. Anal. (2023)] the numerical properties of
this method were explored. These works considered deformations driven
predominantly by bending. Thus, a bending energy with a metric constraint was
considered. We extend these results to the case of an energy with both a
bending component and a nonconvex stretching component, and we also consider
folding across a crease. The proposed discretization of this energy features a
continuous finite element space, as well as a discrete Hessian operator. We
establish the $\Gamma$-convergence of the discrete to the continuous energy and
also present an energy-decreasing gradient flow for finding critical points of
the discrete energy. Finally, we provide numerical simulations illustrating the
convergence of minimizers and the capabilities of the model.
