Computational cardiology is a complex discipline, and one of the cornerstones of every mathematical model of the heart is its fibers. These are usually computed with inexpensive computational models, known as Laplace Dirichlet Rule Based Models (LD-RBM), which yield a very good match with clinical observations. In this work, we obtain a unique system for cardiac fibers which does not rely on potentials, and thus yields what could be considered as the biventricular fiber PDE. The obtained equations are a particular case of the well-known Oseen-Frank equations from liquid crystals theory, which are result in a vector diffusion equation with a unit norm constraint, for which there is a vast literature in the mathematical analysis community. The numerical approximation of these equations is typically done using a Mixed Finite Elements formulation, which has proved ineffective for our application. We have instead devised a projection gradient descent scheme that is robust and scalable, which makes it perfect for High Performance Computing scenarios with high resolution images of the heart. In addition, our approach does not require the computation of function gradients, meaning that it allows us to compute fiber fields of arbitrary accuracy by means of higher order methods, making it a preferable strategy for thin volumes, such as the right ventricle and the atria.