# Area of 8-node rectangular serendipity finite element

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));

• It seems to me that you are trying to use an isoparametric mapping here where you are mapping from the reference cell to the real cell using the eight shape functions of the serendipity element. That seems reasonable because it should map edges to straight lines. But why do you use 8 quadrature points? Don't you just want to use $2\times 2$ Gauss quadrature? Commented Sep 2, 2021 at 1:45
• On second thought, the serendipity element has quadratic functions, so you need a $3\times 3$ Gauss formula and that will then have 9 quadrature points. In any case, however, the 8 you have don't seem reasonable to me. Commented Sep 2, 2021 at 17:39
• @WolfgangBangerth Thanks for your comments. The analytical integral seems to give the correct result. (See update section of my question) Commented Sep 2, 2021 at 23:28
• Sure, you can do that as well. Commented Sep 3, 2021 at 0:30
• The choice of shape functions and the choice of quadrature formula (for integration) is independent. In your case, the 8 serendipity shape functions are quadratic functions, so you need a $3\times 3$ Gauss formula. You can look up the quadrature points and weights in any number of books. Commented Sep 4, 2021 at 19:30

Where I went wrong was that I did not consider the integration at the points $$\xi = 0$$ and $$\eta = 0$$, because I thought they only applied to 9-node elements.
I get the area 4 square units when I account for the integration point at the center of the element, weighted by $$\frac{64}{81}$$