Contact is one of the most challenging nonlinearities to solve in solid mechanics. In traditional linear finite element analysis, the contact surface is only C^0 continuous, as a result, the normal to the contact surface is not continuous. The normal contact force is directed along the normal in the direction of the contact surface, and therefore, the contact force is discontinuous. This issue is tackled in linear finite element analysis using various surface smoothing techniques, however, a better solution is to use isogeometric analysis where the solution space is spanned by smooth spline basis functions. Unfortunately, spline-based isogeometric contact analysis still has limited applicability to industrial computer aided design (CAD) representations. Building analysis suitable mesh from the industrial CAD representations has been a major bottleneck of the computer aided engineering workflow. One promising alternative field of study, intended to address this challenge, is called the immersed finite element method. In this method, the original CAD domain is immersed in a rectilinear grid called the background mesh. This cuts down the model preparation and the mesh generation time from the original CAD domain, but the method suffers from limited accuracy issues. In this dissertation, the original CAD domain is immersed in an envelope domain which can be of arbitrary topological and geometric complexity and can approximate none, some or all of the features of the original CAD domain. Therefore, the method, called the flex representation method, is much more flexible than the traditional immersed finite element method. Within the framework of the flex representation method, a robust and accurate contact search algorithm is developed, that efficiently computes the collision points between the contacting surfaces in a discrete setting. With this information at hand, a penalty based formulation is derived to enforce the contact constraint weakly for multibody and self-contact problems. In addition, the contact algorithm is used to solve various proof-of-concept academic problems and some real world industrial problems to demonstrate the validity and robustness of the algorithms.