clear; close all
%
%  Test of triangular plate bending element
%  Al plate with one edge clamped, load on corner
%
E = 10.e6;      %% E
nu = 0.3;
rho = 0.1;      %% lb/in^3
th = 0.25;      %% thickness
P = 50;         %% applied load
snodes = [0 0
          25 0
          25 15
          0 15];
nsnodes = size(snodes,1);
edges = [1  2
         2  3
         3  4
         4  1];
holes = [12.5 5 2.2];
nholes = size(holes,1);
for i=1:nholes
    [hsnodes,hedges] = hole(holes(i,1:2),holes(i,3),30,nsnodes);
    snodes = [snodes; hsnodes];
    edges = [edges; hedges];
end
hdata.hmax = 3;
[nodes,elems] = mesh2d(snodes,edges,hdata);
nnodes = size(nodes,1);
nelems = size(elems,1);
nodes(:,3) = zeros(nnodes,1);   %% z coord = 0;
figure
for i=1:nelems
    nn = elems(i,[1:3 1]);
    plot3(nodes(nn,1),nodes(nn,2),nodes(nn,3));
    hold on
end
axis equal
set(gca,'XTick',[],'YTick',[],'ZTick',[])
bnodes = findedge(nodes(:,1:2),snodes,[3 4]);
bnodes = find(bnodes==1);
nbnodes = length(bnodes);
bc = zeros(0,2);
for i=1:nbnodes
    for j=1:5  %% constrain all but rotn about z
        bc = [bc; bnodes(i) j];
    end
end
for i=1:nnodes
    bc = [bc; i 6];  %% rotn about z (dof #6) at all nodes
end
bcdof = 6*bc(:,1) -6 + bc(:,2);
ndof = 6*nnodes;
freedof = 1:ndof;
freedof(bcdof) = [];
%
%  Generate stiffness
%
%%bigk = zeros(ndof,ndof);
fprintf('Generating stiffness ... ');
bigk = spalloc(ndof,ndof,5*ndof);
wt = 0;
for i=1:nelems
    nn = elems(i,1:3);
    [k,m,f,dof] = quad4(nn,nodes,E,nu,th,rho);
    bigk(dof,dof) = bigk(dof,dof) + k;
    uuu = zeros(18,1);
    uuu([3 9 15]) = 1;
    wt = wt + sum(m*uuu);
end
fprintf('done\n');
bigd = zeros(ndof,1);
side12 = findedge(nodes(:,1:2),snodes,[1 2]); side12 = find(side12>0);
side41 = findedge(nodes(:,1:2),snodes,[4 1]); side41 = find(side41>0);
corner = intersect(side12,side41);
load = [corner 3 -1000];   %% 1000 lb down at corner
loaddof = 6*corner -6 +3;
bigp = zeros(ndof,1);
bigp(loaddof) = -P;
%arrow(nodes(corner,1:3)+[0 0 .5],nodes(corner,1:3));
fprintf('Solving eqns ... ');
bigd(freedof) = bigk(freedof,freedof)\bigp(freedof);
fprintf('done.\n');
%
%  Plot deformed
%
figure
scale = 10;
for i=1:nelems
    nn = elems(i,[1 2 3 1]);
    xx = nodes(nn,1) + scale*bigd(6*nn-5);
    yy = nodes(nn,2) + scale*bigd(6*nn-4);
    zz = nodes(nn,3) + scale*bigd(6*nn-3);
    plot3(xx,yy,zz);
    patch(xx,yy,zz,-bigd(6*nn-3));  %% w displ = 3rd dof
    hold on
end
dof = [6*corner-5 6*corner-4 6*corner-3];
ha = arrow(nodes(corner,:),nodes(corner,:)+scale*bigd(dof)');
grid
axis equal
set(gca,'XTick',[],'YTick',[],'ZTick',[])
mesh = sprintf('Mesh size %4.1f',hdata.hmax);
deeohf = sprintf('%d free DOF',length(freedof));
result = sprintf('displ. under load %6.3f in',max(-bigd(3:6:ndof)));
title([mesh ', ' deeohf ', ' result])