Adhesive joints/interfaces are increasingly used in assembly of structural components. Most modern adhesives are reinforced at the micron or nm scales to achieve desired mechanical or secondary properties. Unfortunately, adhesive joints are often the weak link in the load bearing capacity of bonded structures. Thus, predicting how the morphology of reinforced adhesives affects the failure response of the bonded structure is critical for design and safety assessment. Direct Numerical Modeling (DNM), which captures all relevant length scales in a single simulation, is one method for accurately predicting the failure response of heterogeneous interfaces at multiple scales. However, DNM is often too computationally expensive, even for modern and next-generation supercomputers. This dissertation presents a high performance numerical framework that uses Computational Homogenization (CH) to model the multi-scale failure response of heterogeneous interfaces in the 3D finite strain setting. The CH formulation assumes a separation of scales, and locally attaches a Representative Unit Cell (RUC) to each point on the macroscopic interface. The different scales are linked through the variational energy equivalence (Hill's lemma) for interfaces, and the macroscopic interface traction-separation response is computationally derived from the micro-scale RUC. In this fashion, the overall macroscopic response is computed with similar accuracy to DNM, but at reduced computational cost. Fully coupled multi-scale simulations that are capable of capturing both the macroscopic response and micro-scale failure processes of bonded structures are important to the "Virtual Materials Testing" paradigm. This dissertation implements the CH formulation in a hierarchically parallel solver that enables simulation of bonded structures that are O(10) mm at the macro-scale, with numerical resolution down to O(10) nm (O(10^6) length scales) within the micro-scale RUCs. Predictive fully coupled simulations containing up to ~54 billion finite elements and over 28 billion nonlinear degrees of freedom are efficiently computed using up to ~400 thousand computing cores. Moreover, large simulations containing up to 1.1 billion finite elements remain possible using the CH solver on only ~1500 computing cores. The excellent scalability of the CH solver, demonstrated using up to 262,144 computing cores, makes it attractive for effectively utilizing next-generation supercomputers as well.Critical to the CH formulation is development of a RUC that accurately describes material morphology and physical response. This dissertation presents an image-based reconstruction procedure for heterogeneous layers that preserves the higher-order statistical description of a material sample. Micro-scale simulations show that reconstructed RUCs that resolve the material statistics optimally balance accuracy of the physical response and size of the computational domain.Quantitatively describing the micro-scale failure process is important to developing better understanding of how microstructure effects the macro-scale behavior of bonded systems. Two metrics are proposed that describe the extent and character of micro-scale failure and its evolution. The metrics are able to help explain complex trends, such as a non-monotonic particle size dependence of the macro-scale fracture toughness. Furthermore, they can illuminate aspects of the complex coupling between macro-scale phenomena and micro-scale failure.