  We present an algorithm in which the all-order strong coupling expansion of
the Abelian U(1) gauge theory with Wilson plaquette action is sampled. In
addition to the vacuum closed surface graphs of the partition function we
propose to also allow for a class of defects (boundaries) related to Wilson
loops in the ensemble. The efficiency of our scheme in estimating various
observables is compared to a standard Metropolis algorithm.
