A primitive variable mixed order formulation of finite element method for solving three dimensional incompressible Navier-Stokes equations is presented. The method of weighted residuals is used for obtaining the approximate solutions of linear and nonlinear partial differential equations. The Physical domain is discretized by using unstructured tetrahedral elements. Unequal order interpolation functions are used for pressure & velocity variables while the temporal discretization is carried out by using an implicit time marching scheme based on finite differencing. One of the major diffculties arising during the finite element solution of an incompressible Navier-Stokes equations is the efficient factorization/preconditioning of the resulting indefinite stiffness matrix. In this work, the formation of an indefinite matrix is avoided by using a pseudo compressibility technique in which an artificial term is introduced into the mass matrix. The artificial term is time dependent and disposed at a later stage once the steady state is reached. Using this approach, the resulting system of equations can then be solved iteratively with standard preconditioners. The non-linear convective term in the Navier-Stokes equations is linearized in time. To diffuse the numerical oscillations which may occur in convection dominated flows, second-orderTaylor-Galerkinstabilization technique is used. The entire solution procedure is encoded in C++ using object oriented programming. One of the special features of this FEM code is that it uses the exact integrals of the shape functions in order to improve the accuracy of the solution, as supposed to any numerical integration schemes. The solution procedure is validated using the benchmark computations for 3D steady incompressible flows.