function [ke, me] = plate_fem_element(dx, dy, params, P)
% PLATE_FEM_ELEMENT  Stiffness and mass of one Q4 Mindlin-Reissner element.
%
%   [ke, me] = plate_fem_element(dx, dy, params, P) integrates, with 2 x 2
%   Gauss points, the bending energy of the rotations and the transverse-
%   shear energy of a rectangular element dx x dy with three unknowns per
%   node (w, theta_x, theta_y):
%
%     U = 1/2 integral ( kappa' Db kappa + gamma' Ds gamma ) dA,
%     kappa = (d theta_x / dx, d theta_y / dy, d theta_x / dy + d theta_y / dx),
%     gamma = (w_x + theta_x, w_y + theta_y),
%     Db = D [1 nu 0; nu 1 0; 0 0 (1 - nu) / 2],   Ds = k G h I,  k = 5/6,
%
%   and the kinetic energy with the translational mass rho h and the rotary
%   inertia rho h^3 / 12.  Full integration of the shear term is what the
%   simulator this lab follows does; on a thin plate it locks, which is why
%   that simulator needs its 80 x 80 mesh to converge (see the README).
%   Degrees of freedom are ordered (w, theta_x, theta_y) per node, nodes
%   anticlockwise from the lower-left corner.

mat = params.material;
h = params.plate.thickness;
nu = mat.nu;
bend = P.D * [1, nu, 0; nu, 1, 0; 0, 0, 0.5 * (1 - nu)];
shear_modulus = mat.E / (2 * (1 + nu));
kappa = 5 / 6;
shear = kappa * shear_modulus * h * eye(2);
rho = mat.rho;
ke = zeros(12, 12);
me = zeros(12, 12);
gp = 1 / sqrt(3);
detj = dx * dy / 4;
for xi = [-gp, gp]
    for eta = [-gp, gp]
        [n, dxi, deta] = plate_fem_shape(xi, eta);
        dnx = dxi * 2 / dx;
        dny = deta * 2 / dy;
        bb = zeros(3, 12);
        bs = zeros(2, 12);
        for a = 1:4
            w = 3 * a - 2; rx = 3 * a - 1; ry = 3 * a;
            bb(1, rx) = dnx(a);
            bb(2, ry) = dny(a);
            bb(3, rx) = dny(a);
            bb(3, ry) = dnx(a);
            bs(1, w) = dnx(a); bs(1, rx) = n(a);
            bs(2, w) = dny(a); bs(2, ry) = n(a);
        end
        ke = ke + (bb' * bend * bb + bs' * shear * bs) * detj;
        for a = 1:4
            for b = 1:4
                nn = n(a) * n(b) * detj;
                me(3 * a - 2, 3 * b - 2) = me(3 * a - 2, 3 * b - 2) + rho * h * nn;
                rotary = rho * h^3 / 12 * nn;
                me(3 * a - 1, 3 * b - 1) = me(3 * a - 1, 3 * b - 1) + rotary;
                me(3 * a, 3 * b) = me(3 * a, 3 * b) + rotary;
            end
        end
    end
end
end