Solve linear elasticity problem in Matlab

dirichlet_example.m

close all;
clear all;

% define geometry
g=[ 2.0000    2.0000    2.0000    2.0000    2.0000    2.0000    2.0000    2.0000
   -0.7169   -0.1488    0.0706    0.6718    0.5185    0.2480   -0.1849   -0.3201
   -0.1488    0.0706    0.6718    0.5185    0.2480   -0.1849   -0.3201   -0.7169
   -0.0316    0.4644    0.0376    0.5215   -0.6358   -0.2630   -0.7560   -0.0676
    0.4644    0.0376    0.5215   -0.6358   -0.2630   -0.7560   -0.0676   -0.0316
         0         0         0         0         0         0         0         0
    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000];

% generate mesh
[p,e,t]=initmesh(g,'Hmax',0.05);

% plot mesh
figure;
line([p(1,t(1,:));p(1,t(2,:))],[p(2,t(1,:));p(2,t(2,:))]);
hold on;
line([p(1,t(2,:));p(1,t(3,:))],[p(2,t(2,:));p(2,t(3,:))]);
line([p(1,t(1,:));p(1,t(3,:))],[p(2,t(1,:));p(2,t(3,:))]);

% initialize matrices
nv=size(p,2);
nt=size(t,2);

K=sparse(nv,nv);
M=sparse(nv,nv);

% quadrature rule
qp=[0.3333    0.2000    0.6000    0.2000
    0.3333    0.6000    0.2000    0.2000];
qw=[-0.2813
    0.2604
    0.2604
    0.2604];

% basis functions
phi{1}=1-qp(1,:)-qp(2,:);
phi{2}=qp(1,:);
phi{3}=qp(2,:);

% basis function gradients
gradhat_phi{1}=repmat([-1;-1],1,length(qw));
gradhat_phi{2}=repmat([1;0],1,length(qw));
gradhat_phi{3}=repmat([0;1],1,length(qw));

% affine mappings
A=cell(2,2);
A{1,1}=p(1,t(2,:))'-p(1,t(1,:))';
A{1,2}=p(1,t(3,:))'-p(1,t(1,:))';
A{2,1}=p(2,t(2,:))'-p(2,t(1,:))';
A{2,2}=p(2,t(3,:))'-p(2,t(1,:))';

% determinants of affine mappings
detA=A{1,1}.*A{2,2}-A{1,2}.*A{2,1};

% inverse mappings
invAt=cell(2,2);
invAt{1,1}=A{2,2}./detA;
invAt{2,1}=-A{1,2}./detA;
invAt{1,2}=-A{2,1}./detA;
invAt{2,2}=A{1,1}./detA;

for i=1:3
    % function values
    u=repmat(phi{i},nt,1);
    ux=bsxfun(@times,invAt{1,1},gradhat_phi{i}(1,:))+bsxfun(@times,invAt{1,2},gradhat_phi{i}(2,:));
    uy=bsxfun(@times,invAt{2,1},gradhat_phi{i}(1,:))+bsxfun(@times,invAt{2,2},gradhat_phi{i}(2,:));

    for j=1:3
        v=repmat(phi{j},nt,1);
        vx=bsxfun(@times,invAt{1,1},gradhat_phi{j}(1,:))+bsxfun(@times,invAt{1,2},gradhat_phi{j}(2,:));
        vy=bsxfun(@times,invAt{2,1},gradhat_phi{j}(1,:))+bsxfun(@times,invAt{2,2},gradhat_phi{j}(2,:));

        K=K+sparse(t(j,:),t(i,:),((ux.*vx+uy.*vy)*qw).*abs(detA),nv,nv);
        M=M+sparse(t(j,:),t(i,:),((u.*v)*qw).*abs(detA),nv,nv);
    end
end

% boundary and interior nodes
D=sort(e(1,:));
I=setdiff(1:nv,D);

% solve eigenvectors
[U,L]=eigs(K(I,I),M(I,I),5,'sm');

% plot eigenvectors
for itr=5:-1:1
    u=zeros(nv,1);
    u(I)=U(:,itr);
    figure;
    subplot(2,1,1);
    pdecont(p,t,u,[0 0]);
    subplot(2,1,2);
    pdesurf(p,t,u);
end