Atrial fibrillation (AF) is the most common arrhythmia encountered clinically, and as the population ages, its prevalence is increasing. Although the CHA2DS2− VASc score is the most used risk-stratification system for stroke risk in AF, it lacks personalization. Patient-specific computer models of the atria can facilitate personalized risk assessment and treatment planning. However, a challenge faced in creating such models is the complexity of the atrial muscle arrangement and its influence on the atrial fiber architecture. This work proposes a semi-automated rule-based algorithm to generate the local fiber orientation in the left atrium (LA). We use the solutions of several harmonic equations to decompose the LA anatomy into subregions. Solution gradients define a two-layer fiber field in each subregion. The robustness of our approach is demonstrated by recreating the fiber orientation on nine models of the LA obtained from AF patients who underwent WATCHMAN device implantation. This cohort of patients encompasses a variety of morphology variants of the left atrium, both in terms of the left atrial appendages (LAAs) and the number of pulmonary veins (PVs). We test the fiber construction algorithm by performing electrophysiology (EP) simulations. Furthermore, this study is the first to compare its results with other rule-based algorithms for the LA fiber architecture definition available in the literature. This analysis suggests that a multi-layer fiber architecture is important to capture complex electrical activation patterns. A notable advantage of our approach is the ability to reconstruct the main LA fiber bundles in a variety of morphologies while solving for a small number of harmonic fields, leading to a comparatively straightforward and reproducible approach.