Fluid-structure interaction (FSI) has gained significant traction in recent decades, with applications spanning various disciplines, including geophysics and biomedicine. In FSI, computational techniques are defined by the choice of discrete domain representation, falling into two main categories: ``boundary-fitted'' or ``non-boundary-fitted'' meshes. Boundary-fitted methods offer high accuracy, but their viability is limited in the presence of large solid deformations. Non-boundary-fitted methods maintain separate and non-matching fluid and structure meshes. Here, structure and fluid are described within a Lagrangian and Eulerian framework, respectively. However, higher mesh resolution is to maintain comparable accuracy, making parallel computing necessary. We present an immersed domain approach for the numerical solution of fluid-structure-contact-interaction (FSCI) problems. The fluid and structure are coupled within the whole overlapping volume, while the different structures in contact are coupled on their surfaces. These couplings are achieved with the method of dual Lagrange multipliers. The nonlinear solution procedure is achieved by solving a sequence of the statically condensed system where only the fluid variables are unknown. We show our general algorithmic framework with our primary parallel computing tools and present a fully coupled simulations for biomedical and industrial applications.