This thesis develops novel methods that combine finite element methods (FEMs) and neural networks (NNs) to solve complex problems constrained by partial differential equations (PDEs). Compared to standard NN solvers, this approach improves accuracy and stability while addressing challenges like boundary conditions and complex geometries. The core idea is interpolating NNs onto FEM spaces, avoiding costly derivatives and ensuring exact integration. Adaptive mesh refinement is introduced for sharper solutions, and the method extends to various function spaces for wider applications. Tested on numerous PDEs, the framework outperforms FEMs in accuracy and offering reliable solutions for inverse problems.