  Many important physical problems, such as fluid structure interaction or
conjugate heat transfer, require numerical methods that compute boundary
derivatives or fluxes to high accuracy. This paper proposes a novel alternative
to calculating accurate approximations of boundary derivatives of elliptic
problems: instead of postprocessing, we describe a new continuous finite
element method based on p-refinement of cells adjacent to the boundary to
increase the approximation order of the derivative on the boundary itself. We
prove that the order of the approximation on the p-refined cells is, in 1D,
determined by the rate of convergence at the knot connecting the higher and
lower order cells and that this idea can be extended, in some simple settings,
to 2D problems. We verify this rate of convergence numerically with a series of
experiments in both 1D and 2D.
