%% clean up
clear
close all

%% load data
exampleCase = 'block';
% Load mesh
fid=fopen([exampleCase,'.grid'],'r');
num=fscanf(fid,'%i',2);
np=num(1);
nt=num(2);
ptmp=fscanf(fid,'%f',[2,np]);
ttmp=fscanf(fid,'%f',[3,nt]);
fclose(fid);
% Change for Matlab
nodes=ptmp';
triangles=ttmp';
clear ptmp ttmp
% Load Dirichlet boundary conditions
DirichletNodes = load([exampleCase,'.DirichletBC']);
% Load Neumann boundary conditions
NeumannSegments = [];
if exist([exampleCase,'.nonzeroNeumannBC'],'file')
    NeumannSegments = load([exampleCase,'.nonzeroNeumannBC']);
end
% Load volume forces
f = zeros(2*np,1);
if exist([exampleCase,'.f'],'file')
    f = [zeros(np,1);-load([exampleCase,'.f'])];
end
% Load elastic parameters
params = ones(size(triangles,1),2);
if exist([exampleCase,'.lm'],'file')
    params = load([exampleCase,'.lm']);
end
% Set lambda and mu from file:
lambda = params(:,1);
mu = params(:,2);
% or set manually...

%% plot
% Plot mesh
triplot(triangles,nodes(:,1),nodes(:,2));
axis equal;
hold on;
% Plot Dirichlet nodes
plot(nodes(DirichletNodes,1),nodes(DirichletNodes,2),'r.');
% Plot Neumann boundary conditions
for i = 1:size(NeumannSegments,1)
    center = (nodes(NeumannSegments(i,1),:)+nodes(NeumannSegments(i,2),:))/2;
    arrowEnd = center + .1*NeumannSegments(i,3:4);
    coords = [center;arrowEnd]';
    axis equal;
    plot(coords(1,:),coords(2,:),'k');
end
