Structural health monitoring (SHM) is of significant importance in the operation of engineering systems to ensure the durability and reliability. In this article, we introduce a Bayesian optimization method using a multioutput Gaussian process to solve the structural fault diagnosis problem. This method utilizes a high fidelity finite element model (FE) of the structure and the impedance/admittance measurements from the structure to identify the location and severity of the damage. The method improves the accuracy of the damage diagnosis by adopting a multioutput Gaussian process as the surrogate model for the full FE model and Thompson sampling approach is used to guide the search for the structural damage in the Bayesian optimization. The detailed algorithms are presented, and the convergence analysis of the method is conducted. We apply our proposed method on simulated synthetic functions and it achieves better performance and higher convergence speed than the traditional mixed input optimization methods. We then apply our method on a real world structural damage identification problem using measured piezoelectric admittance data and illustrate the effectiveness of the proposed method.