We study the space of \(C^{2}\)-smooth isogeometric functions on bilinearly parameterized multi-patch domains \(\Omega \subset \mathbb{R}^{2}\), where the graph of each isogeometric function is a multi-patch spline surface of bidegree \((d,d)\), \(d \in \{5,6 \}\). The space is fully characterized by the equivalence of the \(C^2\)-smoothness of an isogeometric function and the \(G^2\)-smoothness of its graph surface, cf. (Groisser and Peters,2015; Kapl et al.,2015). This is the reason to call its functions \(C^{2}\)-smooth geometrically continuous isogeometric functions. In particular, the dimension of this \(C^{2}\)-smooth isogeometric space is investigated. The study is based on the decomposition of the space into three subspaces and is an extension of the work (Kapl and Vitrih, 2017) to the multi-patch case. In addition, we present an algorithm for the construction of a basis, and use the resulting globally \(C^{2}\)-smooth functions for numerical experiments, such as performing \(L^{2}\) approximation and solving triharmonic equation, on bilinear multi-patch domains. The numerical results indicate optimal approximation order.