Operator networks have emerged as promising deep learning tools for approximating the solution to partial differential equations (PDEs). These networks map input functions that describe material properties, forcing functions and boundary data to the solution of a PDE, i.e., they learn the solution operator of the PDE. In this talk, we consider a new type of operator network called VarMiON, that mimics the variational or weak formulation of PDEs. A precise error analysis of the VarMiON solution reveals that the approximation error contains contributions from the error in the training data, the training error, quadrature error in sampling input and output functions, and a “covering error’’ that measures how well the input training functions cover the space of input functions. Further, the error estimate also depends on the stability constants for the true solution operator and its VarMiON approximation. Numerical experiments are presented for a canonical elliptic PDE to demonstrate the efficacy and robustness of the VarMiON.