This paper presents an explicit finite element integration scheme to compute the stiffness matrices for linear convex quadrilaterals. Finite element formulationals express stiffness matrices as double integrals of the products of global derivatives. These integrals can be shown to depend on triple products of the geometric properties matrix and the matrix of integrals containing the rational functions with polynomial numerators and linear denominator in bivariates as integrands over a 2-square. These integrals are computed explicitely by using symbolic mathematics capabilities of MATLAB. The proposed explicit finite element integration scheme can be applied to solve boundary value problems in continuum mechanics over convex polygonal domains.We have also developed an automatic all quadrilateral mesh generation technique for convex polygonal domain which provides the nodal coordinates and element connectivity.We have demonstrated the proposed explicit integration scheme to solve the Poisson Boundary Value Problem for a linear elastic torsion of a non-circular bar with cross sections having profiles of equilateral triangle, a square and regular polygons (pentagon(5-gon)to icosagon(20gon)) which are inscribed in a circle of unit radius. Monotonic convergence from below is observed with known analytical solutions for the Prandtl stress function and torsional constant .We have shown the solutions in Tables which list both the FEM and exact solutions. The graphical solutions of contour level curves and the corresponding finite element meshes are also displayed.