Aiming to construct a less intrusive contact solver, we propose an accelerated \emph{flexibility method} [1]. This approach constructs a Steklov-Poincaré operator Sc relating contact displacements to contact pressures, reformulating the Signorini-Moreau contact problem as an auxiliary Linear Complementarity Problem (LCP) solved by dedicated algorithms. The original method, which relies on this fully populated operator, has two bottlenecks. First, the construction of this matrix is computationally intensive. Second, it requires substantial storage memory, and the manipulation of such a dense operator is also computationally expensive. Therefore, to revive this method, we explore efficient construction strategies for Sc using both finite element (FEM) and boundary element (BEM) methods. For FEM, we compare Schur complement condensation versus an adapted direct sampling approach. For BEM, we derive Sc from the discretized Somigliana identity, reducing problem dimensionality (3D to 2D). To overcome the difficulties associated with the fully populated operator Sc, we employ hierarchical matrix (H-matrix) [2] approximations via Adaptive Cross Approximation (ACA), reducing computational complexity from O(N^2) to O(N log(N)) and considerably reducing memory requirements. The focus is on deriving a non-intrusive yet efficient computational framework for rigid-to-deformable, deformable-to-deformable, as well as frictionless and frictional contact problems. This framework preserves existing FEM architectures while enabling scalable 3D contact simulations, especially when quasi-static contact problems are solved under varying boundary conditions. Future developments include deformable-deformable contact, friction, and interfacing with open-source codes (FEniCSx, MFEM, in-house A-set and others). [1] Francavilla A., Zienkiewicz O.C., A note on numerical computation of elastic contact problems, Int. J. Numer. Methods Eng., Vol. 9, pp. 913–924, 1975. [2] Bebendorf M., Hierarchical Matrices, Springer, 2008.