In this 3D fracture-based topology optimization framework, the author extended the 2D framework on tailoring fracture resistance for brittle materials (Zhang et al., 2022) to 3D. In the optimization, the topology is described by radial basis functions interpolated level set function, and the problem is solved by the Interface-enriched Generalized Finite Element Method (IGFEM). Cracks are assumed to exist on enriched nodes that are added on the boundary of the geometry to increase the accuracy of approximation. The first part of the work assumes cracks to be semi-circular with the crack plane perpendicular to the boundary, and the crack opening direction to be either parallel to the 𝑋𝑍- or 𝑌𝑍-plane of the global coordinates. An extended framework that assumes crack opening direction perpendicular to the surface first principal stress is also developed at the end. The energy release rates (ERRs) of the cracks are evaluated with the topological derivative method, which requires only a stress field of the geometry and weight functions that relate the stress and stress intensity factors (SIFs). The weight functions are found by a finite element analysis on a cuboid with a crack. This approach is computationally efficient because it eliminates the need of actually modeling and meshing the crack planes in the geometry during optimization. Moreover, a 3D stress recovery technique (or stress improvement procedure, SIP) is used to recover the nodal stress non-locally to improve the accuracy. The objective function is then established with the 𝑝-mean aggregation of the ERRs. Finally, Numerical examples in 3D, including the famous L-bracket benchmark problem, are performed to prove the correctness and capacity of the framework. In conclusion, this extended framework shows more flexibility and provides more information to the optimized design than the 2D framework by considering an added dimension both in analyzed geometry and crack shape, i.e., the effect of the anisotropy of cracks can be captured.