function PE = ComputePE(p)
PE = zeros(length(p.x),1);
%for j = 1:p.dof
j = 3;
  PE = PE - (p.Forces{1}(:,3)+p.Forces{2}(:,3)).*p.x(:,j)/2;
%endplot
