Branched discontinuities appear in problems like crack propagation in brittle materials, stress corrosion cracks, high velocity impact and modeling of polycrystalline grains. The analysis of this class of problems using the classical finite element method (FEM) encounters several difficulties. The singularities at crack fronts and triple joints of polycrystalline grains require strongly refined finite element meshes that must fit the discontinuity surfaces while keeping the aspect ratio of the elements within acceptable bounds. These requirements are not always feasible, specially in three-dimensional simulations, or may require a large number of finite elements. The authors have recently proposed a generalized FEM (GFEM) able to model arbitrary polycrystalline structures  using only the grain topology and a background GFEM discretization that can accommodate grain boundaries and joints arbitrarily located in the mesh. This paper presents extensions of the generalized FEM presented in  that is able to handle arbitrary three-dimensional branched discontinuities and polycrystalline grains while delivering high convergence rates. In the GFEM presented in , the mesh can be refined independently of the discontinuity topology. This, combined with the proposed high-order GFEM approximations, provides a very flexible and robust method that can deliver highly accurate solutions. Some of the capabilities of the method are demonstrated in the solution of three-dimensional branched cracks and in a study of the effect of grain morphology on anelasticity of polycrystals. A series of simulations are carried out on a wide range of grain morphologies in several arrangements of grains. Theoretical analysis and numerical experiments show that the method delivers optimal convergence rates.