This work presents a comprehensive extension of a global–local adaptive mesh strategy for large-scale three-dimensional phase-field fracture simulations. The proposed approach is built upon a goal-oriented error estimator, from which a normalized nodal density field is constructed to guide mesh adaptivity. This formulation enables accurate resolution of highly localized fracture features while maintaining stable mesh quality and significantly reducing computational cost. To efficiently identify critical fracture regions with stringent resolution requirements, a novel global crack sampling strategy based on spatial filtering and clustering is introduced. Subsequently, within Voronoi-partitioned subdomains, a hexagonal close-packing node insertion scheme is employed to achieve efficient node enrichment and smooth, well-structured mesh refinement along complex and evolving crack paths. The final adaptive mesh is generated using a constrained Delaunay tetrahedralization procedure, ensuring conformity to both geometric boundaries and prescribed resolution criteria. To further enhance adaptivity, a self-correcting framework is developed, in which the mapping between estimated error measures and target element sizes is dynamically adjusted at each load increment, thereby providing consistent control of local mesh resolution throughout the simulation. Numerical examples focusing on complex three-dimensional mixed-mode fracture problems demonstrate the robustness, accuracy, and computational efficiency of the proposed method. The results further indicate the potential of the approach for high-fidelity phase-field modeling of complex fracture behavior in large-scale engineering structures.