The stiffness matrix and the nodal forces associated with distributed loads are obtained for a nonhomogeneous anisotropic elastic beam element by the use of complementary energy. The element flexibility matrix is obtained by integrating the complementary-energy density corresponding to six beam equilibrium states, and then inverted and expanded to provide the element-stiffness matrix. Distributed element loads are represented via corresponding internal-force distributions in local equilibrium with the loads. The element formulation does not depend on assumed shape functions and can, in principle, include any variation of cross-sectional properties and load variation, provided that these are integrated with sufficient accuracy in the process. The ability to represent variable cross-sectional properties, coupling from anisotropic materials, and distributed element loads is illustrated by numerical examples.