Modeling of uncertainty within computational mathematics is an active field with many diverse applications. We consider solving partial differential equations with finite element method where some of the data of the problem is uncertain. We explore modeling of uncertain but piecewise constant coefficients of elliptic boundary value problems. The domain of the PDE is bounded and is decomposed into subdomains dependent on the values of coefficients. Either the value or the interface between the subdomains is uncertain. We are interested in understanding the influence of the uncertainty of the coefficients upon the fluxes through the boundaries. We present the model problem and the approaches we used.