function [Area] = Area_Q4_or_Q8 (n_n)   % Revised 2/1/19
% Area of a Q4 or Q8 element by numerical integration

% addpath /clear/www/htdocs/mech517/Akin_FEA_Lib
n_e = 1   ;% number of elements in mesh (for plot)
n_m = n_n ;% number of nodes in mesh (for plot)
n_t = 1   ;% number of element types in mesh (for plot)

if ( n_n == 4 )        ; % number of nodes per element
  x = [0  48.  0  -60.]; % x-coordinates
  y = [0. 64. 100. 40.]; % y-coordinates
  nodes = [1  2  3  4] ; % element connections
  n_q = 1                % number of quadrature pts
  plot_input_2d_mesh (n_e, n_m, n_n, n_t, x, y, nodes);

else                                        ;% Q8 only
  x = [0  48.  0  -60.  31.  29.  -30. -30.];% x-coord
  y = [0. 64. 100. 40.  30.  86.   77.  31.];% y-coord
  %B x = [0  48.  0  -60.  24.  24.  -30. -30.];% straight
  %B y = [0. 64. 100. 40.  32.  82.   70.  20.];% straight
  nodes = [1  2  3  4  5  6  7  8];  % connections
  n_q = 4                            % number of pts
  q8_mesh_plot (x, y, nodes)      ;  % plot curved
end ; % if straight or curved element

% element coordinates as rectangular array
Coord = [x; y]' ; % physical coordinates of nodes

% get the quadrature data over (-1,-1) to (1,1)
[a_q, b_q, w_q] = qp_rule_nat_quad (n_q)  ;

Area = 0.0     ; % initialize area sum
for q = 1: n_q ; % loop over integration points --> --> 
  a = a_q(q) ; b = b_q(q) ; w = w_q(q) ; % get pts & wt

  % Interpolate over (-1,-1) to (1,1)
  [H, DLH] = Serendipity_quads (n_n, a, b) ;
  % H   = element interpolation functions at a,b
  % DLH = interpolation local derivatives at a,b

% calculate the Jacobian, for this point only
  Jacobian = DLH * Coord    ; % numerical 2 by 2 product
  Det_J    = det (Jacobian) ; % numerical determinant
  Area = Area + Det_J * w   ; % increment the area
end ; % loop over integration points <-- <-- <-- <-- <--
% end Area_Q4_or_Q8 =======================================