Thermodynamically consistent phase field theory for multivariant martensitic transformations, which includes large strains and interface stresses, is developed. Theory is formulated in a way that some geometrically nonlinear terms do not disappear in the geometrically linear limit, which in particular allowed us to introduce the expression for the interface stresses consistent with the sharp interface approach. Namely, for the propagating nonequilibrium interface, a structural part of the interface Cauchy stresses reduces to a biaxial tension with the magnitude equal to the temperature-dependent interface energy. Additional elastic and viscous contributions to the interface stresses do not require separate constitutive equations and are determined by solution of the coupled system of phase field and mechanics equations. Ginzburg-Landau equations are derived for the evolution of the order parameters and temperature evolution equation. Boundary conditions for the order parameters include variation of the surface energy during phase transformation. Because elastic energy is defined per unit volume of unloaded (intermediate) configuration, additional contributions to the Ginzburg-Landau equations and the expression for entropy appear, which are important even for small strains. A complete system of equations for fifth- and sixth-degree polynomials in terms of the order parameters is presented in the reference and actual configurations. An analytical solution for the propagating interface and critical martensitic nucleus which includes distribution of components of interface stresses has been found for the sixth-degree polynomial. This required resolving a fundamental problem in the interface and surface science: how to define the Gibbsian dividing surface, i.e., the sharp interface equivalent to the finite-width interface. An unexpected, simple solution was found utilizing the principle of static equivalence. In fact, even two equations for determination of the dividing surface follow from the equivalence of the resultant force and zero-moment condition. For the obtained analytical solution for the propagating interface, both conditions determine the same dividing surface, i.e., the theory is noncontradictory. A similar formalism can be developed for the phase field approach to diffusive phase transformations described by the Cahn-Hilliard equation, twinning, dislocations, fracture, and their interaction.