For many years, torsion of arbitrary cross-sections has been a subject of numerous investigations from theoretical and numerical points of view. As it is well known, the resulting boundary value problem (BVP) governing such phenomenon happens to be a pure Neumann BVP and, therefore, its solutions are determined up to a constant. Among a large plethora of finite element method (FEM) techniques that can be used in this context, most of FEM practitioners resolve this uniqueness issue by fixing the candidate solution to a node of the domain. Although such popular and pinpointing technique is widely spread and works well for practical purposes, it does not have a continuous counterpart and therefore its justification remains a matter of debate. Hence, this self-contained work aims to address the modeling of arbitrary heterogeneous and orthotropic cross-sections as well as the theoretical and numerical aspects of their solutions. In particular, we discuss the existence of weak solutions, well-posedness, regularity of solutions, and convergence of Galerkin’s method for different variational settings (with special focus on a regularized variational approach). Moreover, we establish a connection, at a discrete level, between the convergence of solutions of well-posed variational settings and those solutions coming from the usual practice of fixing a datum at a node. Finally, we discuss some numerical aspects of all the FEM discrete formulations proposed here by performing convergence analysis in L2 and H1 norms. The section of numerical results is closed by presenting a series of study cases ranging from a square cross-section composed of two different materials to an isotropic bridge cross-section for which no analytical solution exists.