clear; close all
%
%  Tapered cantilever beam with holes
%
len = 36;   %% in
base = 15;
tip = 5;
thick = 0.2;
shear = 1000;  %% lb/in on right edge
E = 29.5e6;
nu = 0.3;
rho = 0.283;    %% lb/in^3
elem_size = 1;      %% max element size
%
%  Outer supernodes & edges
%
snodes = [0 -base/2
          len -tip/2
          len tip/2
          0 base/2];
edges = [1  2
         2  3
         3  4
         4  1];
%
%  holes (center x, y, rad)
%
holes = [0.2*len  0  3.2
         0.45*len  0  2.3
         0.7*len  0  1.5
         0.9*len  0  1.0];
nholes = size(holes,1);
inc = 4;
for i=1:nholes
    [ss,ee] = hole(holes(i,1:2),holes(i,3));
    ee = ee+inc;
    inc = inc + size(ss,1);
    snodes = [snodes; ss];
    edges = [edges; ee];
end
hdata.hmax = elem_size;
[nodes,elems] = mesh2d(snodes,edges,hdata);
nnodes = size(nodes,1);
nelems = size(elems,1);
ndof = 2*nnodes;
%
%  Find left edge, constrain, color red
%
ind = find(nodes(:,1) < eps);
bc = zeros(0,2);
for i=1:size(ind,1);
    bc = [bc; [ind(i) 1]; [ind(i) 2] ];
end
bcdof = 2*bc(:,1) -2 + bc(:,2);
hold on
plot(nodes(ind,1),nodes(ind,2),'r*','LineWidth',2);
%
%  Find right edge, apply shear load (force per unit len)
%
ind = find(nodes(:,1) >= len-eps);
plot(nodes(ind,1),nodes(ind,2),'g*','LineWidth',2);
bigp = zeros(ndof,1);
for i=2:length(ind)
    elen = (nodes(ind(i),2) - nodes(ind(i-1),2));
    dof = 2*[ind(i) ind(i-1)];
    bigp(dof) = bigp(dof) - shear*elen/2;
end
bigk = zeros(ndof,ndof);
uuu = [1 0 1 0 1 0]';
wt = 0;
for i=1:nelems
    n1 = elems(i,1);
    n2 = elems(i,2);
    n3 = elems(i,3);
    [k,dof,m] = triak(nodes,n1,n2,n3,thick,E,nu,rho);
    bigk(dof,dof) = bigk(dof,dof) + k;
    wt = wt + sum(m*uuu);
end
title(sprintf('Weight %8.3f lb',wt));
freedof = 1:ndof;
freedof(bcdof) = [];
bigd = zeros(ndof,1);
bigd(freedof) = bigk(freedof,freedof)\bigp(freedof);
%
%  Plot deformed, get stresses
%
fprintf('Hit ENTER to continue'); pause; fprintf('\n');
figure
scale = 50;
color = 'k';
xc = zeros(nelems,1);
yc = zeros(nelems,1);
xd = zeros(nnodes,1);
yd = zeros(nnodes,1);
for i=1:nnodes
    x = nodes(i,1);  y = nodes(i,2);
    dofx = 2*i-1;
    xd(i) = x + scale*bigd(dofx);       % deformed position
    dofy = 2*i;
    yd(i) = y + scale*bigd(dofy);
    xc(i) = mean(nodes([n1 n2 n3],1));  % need centroid for later
    yc(i) = mean(nodes([n1 n2 n3],2));
end;
set(gca,'XTick',[],'XLim',[-1 len+1]);
set(gca,'YTick',[]);
axis equal
sig = zeros(3,nelems);
for i=1:nelems
    n1 = elems(i,1);
    n2 = elems(i,2);
    n3 = elems(i,3);
    %loop = [n1 n2 n3 n1];
    %plot(xd(loop),yd(loop),'--');
    dof = [2*n1-1 2*n1 2*n2-1 2*n2 2*n3-1 2*n3];
    sig(:,i) = trias(nodes,n1,n2,n3,thick,E,nu,bigd(dof));
end
svm = zeros(nnodes,1);
%
%  Average von Mises stresses at nodes & plot
%
for i=1:nnodes
    x = nodes(i,1);
    y = nodes(i,2);
    svm(i) = svm_node(i,x,y,elems,xc,yc,sig);
end;
svmmax = max(svm);
for i=1:nelems
    n1 = elems(i,1);
    n2 = elems(i,2);
    n3 = elems(i,3);
    nn = [n1 n2 n3 n1];
    patch(xd(nn),yd(nn),svm(nn)/svmmax);
end;
title(sprintf('von Mises stresses, ksi (max %5.1f)',svmmax/1000));
axis equal
set(gca,'XTick',[],'YTick',[]);