function out=element_quad4(task,node,coef,order)
% K=element_quad4('stiff',nodes,coef)
% 4-noded isoparametric rectangle 2D
% input:
%   task    'stiff'
%   node   node(:,i) are xyz coordinates of node i
%   coef    the coefficient
%   order   polynomial accuracy
% output:
%   K       local stiffness matrix

% reference quad: [-1,1] x [-1,1]
%
%    4 [-1, 1] --- 3 [1,1] 
%        |
%        |
%        |
%    1 [-1,-1] --- 2 [1,-1]
% basis functions
f1= [1 -1
    -1 1]/4;  % node=[-1,-1] (x1-1)*(x2-1)=-x1*x2-x1-x2+1
f2=[-1 +1
    -1 +1]/4; % node=[1,-1] (x1+1)*(x2-1)=x1*x2-x1+x2-1
f3=[1 1
    1 1]/4; % node=[1,1] (x1+1)*(x2+1)=x1*x2+x1+x2+1
f4=[-1  -1
    +1 +1]/4; % node=[-1,1] (x1-1)*(x2+1)=x1*x2+x1-x2-1
nn=4;      % number of nodes
fun={f1,f2,f3,f4};
nf=size(coef,1); % number of fields
ndof=nn*nf*nf; % dofs in element
for i=1:2, % isoparametric mapping
    map{i}=node(i,1)*f1+node(i,2)*f2+node(i,3)*f3+node(i,4)*f4;
end
if ~exist('order','var')
    order=3;
end
switch order
    case 1 % 1 point Gauss quadrature, exact for polynomials order 1
        gn=[0;0];
        gw=4;
    case {2,3} % 2 point quadrature, exact for pols order 3
        % 2x2 point Gauss quadrature
        s=1/sqrt(3);
        gn=[s,s; s,-s; -s,s; -s,-s]';
        gw=ones(1,4);
    otherwise
        error('bad order')
end
switch task
    case 'stiff'
        % compute local stiffness matrix
        K=stiff_pol(coef,fun,map,gn,gw);
        out=reshape(K,[ndof,ndof]);
    otherwise
        error('bad function')
end