clear;

Lx = 2;
Ly = 3;
Nx = 300;
Ny = 300;

mesh = create_slab([Lx, Ly], [Nx Ny]);
type = 24;

xy = mesh.Nodes(:, 2:3);
node_sel = xy(:,1) > Lx / (10*Nx) & xy(:,1) < Lx - Lx/(10*Nx) & ...
           xy(:,2) > Ly / (10*Ny) & xy(:,2) < Ly - Ly/(10*Ny);
mesh.Nodes(node_sel,2:3) = mesh.Nodes(node_sel,2:3) + ...
    (rand(sum(node_sel), 2)-0.5) * Lx / (6*Nx);
       
tic;
dim = floor(type / 10);
n_elems = Nx * Ny;

[w, N, dNx, x] = fe_quad(mesh.Nodes, mesh.Elements, type, 2);

M_vals = block_product(N.', bsxfun(@times, repmat(N, n_elems, 1), w));
K_vals = multi_block_product( ...
    inner_block_transpose(dNx, size(dNx,1)/n_elems, size(dNx,2)), ...
    bsxfun(@times, dNx, kron(w, ones(dim, 1))));

elem_sel = mesh.Elements(:, 2) == type;
n_nod = mod(type, 10);
node_inds = mesh.Elements(elem_sel, 4 + (1 : n_nod));

rows = repmat(node_inds.', 1, n_nod);
cols = kron(node_inds, ones(n_nod, 1));

M = sparse(rows, cols, M_vals);
K = sparse(rows, cols, K_vals);
t_new = toc;

tic;
[M_old, K_old] = model_mk(mesh);
t_old = toc;