function [Area] = Area_T6_curved (n_q);% Revised 3/17/20
%  Area of a curved T6 element by numerical integration
%  using unit coordinates interpolation and quadrature
%      n_q = number of quadrature points required
%  Area = integral dx dy          (= integral 1 dx dy)
%  Area = integral |J(r,s)| dr ds         
%  If straight with uniform nodes then |J| degree = 0
%  If curved or non-uniform spaces |J| degree = 1+1 = 2
%
%  Area = 0 + Sum_(q=1, n_q)  |J(r_q, s_q)| * w_q
%  For degree = 0, 1, 2, 3, 4, 5,  6,  ...
%  need n_q  >= 1, 1, 3, 4, 6, 7, 12,  ...

n_n = 6                    ;% number of nodes per element
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)

% Define vertices then mid-sides, counter-clockwise (why)
nodes = [1   2   3    4    5     6]       ; % connections
x     = [0.  48. 0.   31.  29.  -20.]         ; % x-coord
y     = [0.  64. 100. 30.  86.   50.]         ; % y-coord 
  
% Set element coordinates as rectangular array
Coord = [x; y]'         ; % physical coordinates of nodes

% plot the curved triangle as a data check
% addpath /clear/www/htdocs/mech517/Akin_FEA_Lib
% t6_mesh_plot (x, y, nodes)           ; % plot curved T6 

% get the quadrature data over (0,0) to (1,0) to (0,1)
OK = [1, 3, 4]         ; % the only valid n_q inputs here
if ( ~any (n_q == OK) )    ; % input not on required list
  fprintf ('Non-existant rule %i requested \n', n_q)
  error   ('n_q invalid in function Area_T6_curved')
end                                       ; % if bad data
%  Allocate storage for quadrature data, and extract them
r_q=zeros(n_q,1); s_q=zeros(n_q,1); w_q=zeros(n_q,1);
if ( n_q == 1 ) % exact for constant or linear polynomial
  r_q = [ 1 / 3.] ; s_q = [ 1 / 3.] ; w_q = [ 1 / 2.] ;

elseif ( n_q == 3 )   % exact for 2-rd degree polynomials
 r_q=[6.6666666666666667e-01, 1.6666666666666667e-01, ...
      1.6666666666666667e-01 ] ;
 s_q=[1.6666666666666667e-01, 6.6666666666666667e-01, ...
      1.6666666666666667e-01 ] ;
 w_q=[1.6666666666666666e-01, 1.6666666666666667e-01, ...
      1.6666666666666667e-01 ] ;

elseif ( n_q == 4 )   % exact for 3-rd degree polynomials
  r_q = [ 1.0, 0.6, 0.6, 1.8 ] / 3.  ;
  s_q = [ 1.0, 0.6, 1.8, 0.6 ] / 3.  ;
  w_q = [ -27, 25,  25,  25  ] / 96. ;
end                          ; % if which rule to extract

Area = 0.0                       ; % initialiaze area sum
for q = 1: n_q ; % loop over integration points ---> ---> 
  r = r_q(q) ; s = s_q(q) ; w = w_q(q)   ; % get pts & wt
  % H   = element interpolation functions at r, s
  % DLH = interpolation local derivatives at r, s
    H = [(1-3*r-3*s+2*r*r+4*r*s+2*s*s), (2*r*r - r), ...
         (2*s*s-s), 4*(r-r*r-r*s), 4*r*s, 4*(s-r*s-s*s)];
    DLH = [(-3+4*r+4*s), (-1+4*r), 0,   ...
           (4-8*r-4*s), (4*s), (-4*s) ; ...       % dH/dr
           (-3+4*r+4*s), 0, (-1+4*s),   ...   
           (-4*r), (4*r), (4-4*r-8*s)] ;          % dH/ds

% calculate the Jacobian, for this point only
  Jacobian = DLH * Coord  ; % 2 by 2 Jacobian matrix at q
  Det_J    = det (Jacobian)   ; % scalar determinate at q
  Area = Area + Det_J * w        ; % increment area value
end ; % loop over integration points  <--- <--- <--- <---
fprintf('Quadrilateral area is %10.4e \n', Area);% result
% end Area_T6_curved ====================================