This paper introduces a rigorous energetic framework for modeling crack propagation in isotropic linear elastic media, established upon a discrete variational formulation. Moving beyond traditional local growth criteria that depend on asymptotic crack-tip field extractions, the proposed formulation is derived from the Principle of Minimum Potential Energy and implemented via the Dual Boundary Element Method. The core innovation is the Potential Energy Minimization Algorithm which conceptualizes crack evolution as a discrete shape-optimization sequence. By determining admissible propagation directions through global energetic stationarity, the algorithm circumvents numerical instabilities and mesh sensitivities associated with explicit Stress Intensity Factor evaluation or crack-tip enrichment strategies. A significant advantage of this approach is the ability to preserve global equilibrium at each incremental step while maintaining fidelity in the presence of stress gradients. The robustness of the PEMA is validated through a suite of benchmark problems, encompassing mixed-mode configurations, stress concentrators, and shielding and amplification effects in multi-crack interactions. The numerical results demonstrate that the PEMA maintains consonance with the Maximum Tangential Stress criterion and experimental data and offers superior trajectory stability. By evaluating the energetic functional through boundary quantities, the framework provides a computationally efficient and consistent foundation for structural integrity assessments and fracture simulations.