We develop and evaluate a method for learning solution operators to nonlinear problems governed by partial differential equations (PDEs). The approach is based on a finite element discretization and aims at representing the solution operator by a multilayer perceptron (MLP) that takes problem data variables as input and gives a prediction of the finite element solution as output. The variables will typically correspond to parameters in a parametrization of input data such as boundary conditions, coefficients, and right-hand sides. The output will be an approximation of the corresponding finite element solution, thus enabling support and enhancement by the standard finite element method (FEM) both theoretically and practically. The loss function is most often an energy functional and we formulate efficient parallelizable training algorithms based on assembling the energy locally on each element. For large problems, the learning process can be made more efficient by using only a small fraction of randomly chosen elements in the mesh in each iteration. The approach is evaluated on several relevant test cases, where learning the finite element solution operator turns out to be beneficial, both in its own right but also by combination with standard FEM theory and software.