Liquid crystal polymeric networks (LCNs) are materials where a nematic liquid crystal is coupled with a rubbery material. When actuated with heat or light, the interaction of the liquid crystal with the rubber creates complex shapes. Thin bodies of LCNs are natural candidates for soft robotics applications. Starting from the classical 3D trace energy formula of Bladon, Warner and Terentjev (1994), we derive a 2D membrane energy as the formal asymptotic limit of the 3D energy and characterize the zero energy deformations. The membrane energy lacks convexity properties, which lead to challenges for the design of a sound numerical method. We discretize the problem with conforming piecewise linear finite elements and add a higher order bending energy regularization to address the lack of convexity. We prove that minimizers of the discrete energy converge to zero energy states of the membrane energy in the spirit of Gamma convergence; this includes the presence of creases. We solve the discrete minimization problem via an energy stable gradient flow scheme. We present computations showing the geometric effects that arise from liquid crystal defects as well as computations of nonisometric origami, both within and beyond theory. This work is joint with R.H. Nochetto and S. Yang.