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 = 10;         %% applied load
snodes = [0 0
          15 0
          15 5
          10 5
          10 12
          0 12];
nsnodes = size(snodes,1);
edges = [1  2
         2  3
         3  4
         4  5
         5  6
         6  1];
holes = [5 4 1.2];
holes = [];
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 = 0.8;
[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,[2 3]);
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);
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
bigd = zeros(ndof,1);
side56 = findedge(nodes(:,1:2),snodes,[5 6]); side56 = find(side56>0);
side61 = findedge(nodes(:,1:2),snodes,[6 1]); side61 = find(side61>0);
corner = intersect(side56,side61);
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));
bigd(freedof) = bigk(freedof,freedof)\bigp(freedof);
%
%  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',[])
weight = sprintf('Weight %5.3f lb',wt);
result = sprintf('displacement under load %6.3f in',max(-bigd(3:6:ndof)));
title([weight ', ' result])