
clear
clf

% Simple finite element solver
%
% Acoustic wave equation  u_tt - div grad(u) = f 
% 1-D regular grid

% Basis function
%                 1 at x = x_i
%    phi_i = 		0 else
%						ax between x_i-1 and x_i, a=1/(x_i-x_i-1)
%                -bx between x and x_i+1, b=1/(x_i+1-x_i)


% We need to solve the system
%
% 			b_i A_ij = g_i				(1)
%
% or
%
%      A_ij^T b_i = g_i				(2)
%
% with the solution
%
%				  b_i = (A_ij^T)^-1 g_i			(3)
%


% Basic parameters

nt=1000;
v=1.;
eps=1;


nx=250;
dx=1/(nx-1);
xx=(0:(nx-1))*dx;
%x=(-cos((0:nx-1)*pi/(nx-1))+1)/2;
x=xx;

% Below is a suggestion how to turn the regular FEM into an irregular grid
% by perturbing the length of each element with a random perturbation with
% max amplitude 0.25*dx You can look at the distribution of grid sizes by
% using the hist function -> hist(h)
%pert=.25;
%x(2:nx-1)=x(2:nx-1)+dx*(rand([1 nx-2])-.5)*2*pert;

% h contains the length of the elements using the difference operation
h=diff(x);

dtfd=.5*eps*dx/v;
dtfe=.5*eps*min(h)/v;

dt=min([dtfd,dtfe]);
disp(sprintf(' dtfd = %g, dtfe = dtfe %g',dtfd,dtfe));

% source function

p=20*dt;
t=(1:nt)*dt; t0=3*p;
src=-1/p^2*(t-t0).*exp(-1/p^2*(t-t0).*(t-t0));
src=diff(src);src(nt)=0.;


% Calculation of the rhs of equation (1)

% source term is a spike at x=nx/2*dx; 

% Quelle
f=zeros([nx 1]);f(nx/2)=1.;
% Randwert links  u_1  int{ nabla phi_1  nabla phij }
u1  = 0;   f(1) = 0;
% Randwert links  u_nx int{ nabla phi_nx nabla phij }
unx = 0;   f(nx)= 0;

% assemble matrix Aij

A=zeros(nx);

for i=2:nx-1,
   for j=2:nx-1,
      if i==j,
         A(i,j)=1/h(i-1)+1/h(i);
      elseif i==j+1
         A(i,j)=-1/h(i-1);
      elseif i+1==j
         A(i,j)=-1/h(i);
      else
         A(i,j)=0;
      end
   end
end

% incorporate velocity in to Matrix A

A=A*v^2;

% assemble matrix Mij

M=zeros(nx);

for i=2:nx-1,
   for j=2:nx-1,
      if i==j, 
         M(i,j)=h(i-1)/3+h(i)/3;
      elseif j==i+1
         M(i,j)=h(i)/6;
      elseif j==i-1
         M(i,j)=h(i)/6;
      else
         M(i,j)=0;
      end
   end
end

% Solve 

% M^T c_tt + A^T c = f 

u=zeros([nx 1]);
uold=u;
unew=u;
p=zeros([nx 1]);
pold=p;
pnew=p;


% Time stepping

for i=1:nt,
   
   % FEM
   unew(2:nx-1) = dt^2 * inv(M(2:nx-1,2:nx-1)) * ...
      				(f(2:nx-1)*src(i)-A(2:nx-1,2:nx-1)*...
      				u(2:nx-1)) + 2*u(2:nx-1) - uold(2:nx-1);
   uold=u;
   u=unew;
   
   
   % FD
   
   lp=v^2*(csh(p,-1)'-2*p+csh(p,1)')/dx^2;
   lp(nx/2)=lp(nx/2)+src(i)/dx;
   pnew=2*p-pold+lp*dt^2;
   pold=p;
   p=pnew;
   p(1)=0;
   p(nx)=0;
   
   % Display both
   plot(xx,p,'b-',x,u,'r-')
   title(' FD - blue; FEM - red ')
   drawnow
end

   

