This paper presents an explicit integration scheme to  compute  the stiffness matrix of twelve node and sixteen node linear convex quadrilateral   finite elements of Serendipity and Lagrange families  using an explicit integration scheme and discretisation of polygonal domain by such finite elements using a novel auto mesh generation technique, In finite element analysis, the boundary value problems governed by  second order linear partial  differential equations, the element stiffness matrices are expressed as integrals of the product of global derivatives over the linear convex quadrilateral region. These matrices can be shown to depend on the material properties matrices and the  matrix  of integrals with integrands as rational functions with polynomial numerator and the linear denominator  (4+)  in the  bivariates    over a   2-square   (-1 ) with the  nodes on the boundary and in the interior of this simple domain. The finite elements  up to cubic order have nodes only on the boundary for Serendipity family and the finite elements with boundary as well as some interior nodes belong to the Lagrange  family. The first order element is the bilinear convex quadrilateral finite element which is an exception and it belongs to both the families. We have for the present ,the  cubic order finite elements which havee 12 boundary nodes  at the nodal coordinates {(-1,-1),(1,-1),(1,1),(-1,1),(-1/3,-1), (1/3,-1),(1,-1/3),(1,1/3),(1/3,1),(-1/3,1),(-1,1/3),(-1,-1/3)} and the four  interoior nodal coordinates at the points (-1/3,-1/3),(1/3,-1/3),(1/3,1/3),(-1/3,1/3)}  in  the local parametric space ( In this paper, we have computed the integrals of  local derivative products  with linear denominator (4+)    in  exact  forms  using the symbolic mathematics capabilities of MATLAB. The proposed explicit  finite  element integration scheme can be then applied to solve boundary value problems  in continuum mechanics over convex polygonal domains. We have  also developed  a novel auto mesh  generation  technique of all 12-node and 16-node linear(straight edge) convex quadrilaterals  for a  polygonal  domain  which provides the nodal coordinates and the element connectivity. We have used  the explicit integration  scheme and this  novel auto mesh generation technique to solve the Poisson equation  u ,where u is an  unknown physical variable  and      in with  Dirichlet  boundary conditions over  the  convex  polygonal domain.