An integral equation approach based on bi-cubic Hermite surfaces has been proposed to solve the three-dimensional complex Stokes flow problems. The numerical solutions are obtained by utilizing the boundary collocation method as well as the continuous distribution of Stokeslets, which are the fundamental solutions of the steady Stokes equations. The quadrilateral surface elements are represented with the bi-cubic Hermite functions that allow the continuous variation of the surface normal vectors between neighbouring elements. The unique normal vectors at the vertices are obtained using the MeanWeighted by Sine and Edge Length Reciprocal (MWSELR) approach. The singular integrations are evaluated using the Tanh - Sinh quadrature rule, meanwhile non-singular integrals are evaluated using the Gauss-Legendre quadrature rule. For the case of freely moving/falling rigid bodies, the integral equations are fully coupled with the equation of motion and the resulting algebraic equations are solved iteratively using the preconditioned Krylov subspace methods provided within the PETSc library. The accuracy of the numerical algorithm is initially verified for the three-dimensional unbounded Stokes flow around a sphere. Then the algorithm is applied to the more complex problems such as the sedimentation of spherical particles and helical swimmers. The numerical simulations show parallel scalability for a large number of processors.