  This paper presents a numerical method for div-curl systems with normal
boundary conditions by using a finite element technique known as primal-dual
weak Galerkin (PDWG). The PDWG finite element scheme for the div-curl system
has two prominent features in that it offers not only an accurate and reliable
numerical solution to the div-curl system under the low $H^\alpha$-regularity
($\alpha>0$) assumption for the true solution, but also an effective
approximation of normal harmonic vector fields regardless the topology of the
domain. Results of seven numerical experiments are presented to demonstrate the
performance of the PDWG algorithm, including one example on the computation of
discrete normal harmonic vector fields.
