Broadly applicable analytical algorithms for workspace of serial manipulators with non-unilateral constraints are developed and illustrated. The Jacobian row-rank deficiency method is employed to determine the singularities of these manipulators. There are four types of singularity sets: Type I: position Jacobian singularities; Type II: instantaneous singularities that are due to a generalized joint that is reaching its apex; Type III: domain boundary singularities, which are associated with the initial and final values of the time interval; Type IV: coupled singularities, which are associated with a relative singular Jacobian, where the null space is reduced in one sub-matrix due to either of two occurrences: a Type II or a Type III singularity. All of the singular surfaces are hypersurfaces that extend internally and externally the workspace envelope. Intersecting singular surfaces identifies singular curves that partition singular surfaces into subsurfaces, and a perturbation method is used to identify regions (curve segments/surface patches) of the hypersurfaces that are on the boundary. The formulation is illustrated by implementing it to a spatial 3-degree of freedom (DOF) and a spatial 4-DOF manipulator.