I am trying to compute the area of an 8-node rectangular serendipity finite element from the equation $$ \sum_{i= 1}^8 det \, J(\xi,\eta) \cdot W_i $$ based on Gaussian quadrature, where $$ J(\xi,\eta) = \begin{vmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \eta} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \xi} \end{vmatrix} \\ x(\xi, \eta) = \sum_{i= 1}^8 N_i^{8Q}(\xi, \eta) \cdot x_i \\ y(\xi, \eta) = \sum_{i= 1}^8 N_i^{8Q}(\xi, \eta) \cdot y_i \\ $$
The shape functions $N_i^{8Q}(\xi, \eta)$ and its derivates are from this appendix; $W_i = \frac{5}{9}$ and the Gauss points used are $\pm \sqrt{\frac{3}{5}}$
Unfortunately I do not get the answer that I expect for a square shape whose, nodal points are at its corners and the midpoints of the edges.
Is Gaussian quadrature exact in 2d or it an approximation?
Update
My later attempt was to avoid Gauss Quadrature and analytically integrate the formula
$$ \int_{-1}^1 \int_{-1}^1 J(\xi,\eta) \, d\eta \, d\xi $$
The resulting expression is much shorter than the result of my attempt to use Gaussian Quadrature, and more importantly gives correct results for 4 different quadrilaterals (square, rhombus, parallelogram, and trapezium). Is this approach is valid?
Update 2
I wrote a simple maxima script to calculate the area of an 8-node serendipity square finite element of side 2 centered at the origin. The area of the element ought to be 4 but I get 8. In fact all four quadrilaterals that I tested have twice the expected area. What's more surprising is that although I initially used Gauss points $\pm \sqrt{\frac{1}{3}}$ (based on this article), the area is still 8 if I use a Gauss point $\pm \sqrt{\frac{3}{5}}$ or any other value.
/* Run in (wx)maxima or online at http://maxima.cesga.es/index.php */
/*(paste and click the button labelled "Clic"*/
N[1](xi,eta) := -1/4*(1 - xi)*(1 - eta)*(1 + xi + eta);
N[2](xi,eta) := 1/2*(1 - xi)*(1 + xi)*(1 - eta);
N[3](xi,eta) := -1/4*(1 + xi)*(1 - eta)*(1 - xi + eta);
N[4](xi,eta) := 1/2*(1 + xi)*(1 + eta)*(1 - eta);
N[5](xi,eta) := -1/4*(1 + xi)*(1 + eta)*(1 - xi - eta);
N[6](xi,eta) := 1/2*(1 - xi)*(1 + xi)*(1 + eta);
N[7](xi,eta) := -1/4*(1 - xi)*(1 + eta)*(1 + xi - eta);
N[8](xi,eta) := 1/2*(1 - xi)*(1 + eta)*(1 - eta);
define(DN1Dxi(xi,eta), diff(N[1](xi,eta), xi, 1));
define(DN2Dxi(xi,eta), diff(N[2](xi,eta), xi, 1));
define(DN3Dxi(xi,eta), diff(N[3](xi,eta), xi, 1));
define(DN4Dxi(xi,eta), diff(N[4](xi,eta), xi, 1));
define(DN5Dxi(xi,eta), diff(N[5](xi,eta), xi, 1));
define(DN6Dxi(xi,eta), diff(N[6](xi,eta), xi, 1));
define(DN7Dxi(xi,eta), diff(N[7](xi,eta), xi, 1));
define(DN8Dxi(xi,eta), diff(N[8](xi,eta), xi, 1));
define(DN1Deta(xi,eta), diff(N[1](xi,eta), eta, 1));
define(DN2Deta(xi,eta), diff(N[2](xi,eta), eta, 1));
define(DN3Deta(xi,eta), diff(N[3](xi,eta), eta, 1));
define(DN4Deta(xi,eta), diff(N[4](xi,eta), eta, 1));
define(DN5Deta(xi,eta), diff(N[5](xi,eta), eta, 1));
define(DN6Deta(xi,eta), diff(N[6](xi,eta), eta, 1));
define(DN7Deta(xi,eta), diff(N[7](xi,eta), eta, 1));
define(DN8Deta(xi,eta), diff(N[8](xi,eta), eta, 1));
f : sqrt(1/3);
/* x-coordinates */
x[1] : -1;
x[2] : -1;
x[3] : -1;
x[4] : 0;
x[5] : 1;
x[6] : 1;
x[7] : 1;
x[8] : 0;
/* y-coordinates */
y[1] : 1;
y[2] : 0;
y[3] : -1;
y[4] : -1;
y[5] : -1;
y[6] : 0;
y[7] : 1;
y[8] : 1;
define(DxDxi(xi,eta) , DN1Dxi(xi, eta)*x[1] + DN2Dxi(xi, eta)*x[2] + DN3Dxi(xi, eta)*x[3] + DN4Dxi(xi, eta)*x[4] + DN5Dxi(xi, eta)*x[5] + DN6Dxi(xi, eta)*x[6] + DN7Dxi(xi, eta)*x[7] + DN8Dxi(xi, eta)*x[8]);
define(DxDeta(xi, eta) , DN1Deta(xi, eta)*x[1] + DN2Deta(xi, eta)*x[2] + DN3Deta(xi, eta)*x[3] + DN4Deta(xi, eta)*x[4] + DN5Deta(xi, eta)*x[5] + DN6Deta(xi, eta)*x[6] + DN7Deta(xi, eta)*x[7] + DN8Deta(xi, eta)*x[8]);
define(DyDxi(xi, eta) , DN1Dxi(xi, eta)*y[1] + DN2Dxi(xi, eta)*y[2] + DN3Dxi(xi, eta)*y[3] + DN4Dxi(xi, eta)*y[4] + DN5Dxi(xi, eta)*y[5] + DN6Dxi(xi, eta)*y[6] + DN7Dxi(xi, eta)*y[7] + DN8Dxi(xi, eta)*y[8]);
define(DyDeta(xi, eta) , DN1Deta(xi, eta)*y[1] + DN2Deta(xi, eta)*y[2] + DN3Deta(xi, eta)*y[3] + DN4Deta(xi, eta)*y[4] + DN5Deta(xi, eta)*y[5] + DN6Deta(xi, eta)*y[6] + DN7Deta(xi, eta)*y[7] + DN8Deta(xi, eta)*y[8]);
/* Jacobian */
define(J(xi,eta), DxDxi(xi,eta)*DyDeta(xi,eta) - DxDeta(xi,eta)*DyDxi(xi,eta));
/* Area */
ratsimp(J(-f,f) + J(-f,0) + J(-f,-f) + J(0,-f) + J(f,-f) + J(f,0) + J(f,f) + J(0,f));
