The capability to model the nonlinear magnetohydrodynamic (MHD) evolution of stellarator plasmas is developed by extending the M3D-C1 code to allow non-axisymmetric domain geometry. We introduce a set of logical coordinates, in which the computational domain is axisymmetric, to utilize the existing finite-element framework of M3D-C1. A C1 coordinate mapping connects the logical domain to the non-axisymmetric physical domain, where we use the M3D-C1 extended MHD models essentially without modifications. We present several numerical verifications on the implementation of this approach, including simulations of the heating, destabilization, and equilibration of a stellarator plasma with strongly anisotropic thermal conductivity, and of the relaxation of stellarator equilibria to integrable and non-integrable magnetic field configurations in realistic geometries.