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(20-gon)) 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.