%mesh
%domain
x_min = 0;
y_min = 0;
x_max = 100;
y_max = 100;
h    =  2;
%mesh by distmesh
fd    =  @(p) drectangle(p,x_min,x_max,y_min,y_max);
[nodes,elements] =  distmesh2d(fd,@huniform,h,[x_min,y_min;x_max,y_max],[x_min,...
               y_min;x_min,y_max;x_max,y_min;x_max,y_max]);
%**************************************************************************
% input parameter
dt            = 0.0006; %the time interval(s)
nt            = 101;
v             = 1000;  %the velocity of acoustic wave(m/s).
ppw           = 20;
num           = length(nodes);
num_elem      = length(elements);

%allocate memory
%%sparse
K             = sparse(num,num);
M             = sparse(num,num);
F             = sparse(num,1);
RES           = zeros(num,nt);
P             = sparse(num,1);
Pold          = sparse(num,1);
Pnew          = sparse(num,1);
II            = speye(num);
distance      = zeros(num,1);
grid_distance = zeros(num_elem,3);

% search the load location,return node number
loca_x = x_max/2;
loca_y = y_max/2;
for ii =1:1:num
    distance(ii,1) = (nodes(ii,1)-loca_x)^2+(nodes(ii,2)-loca_y)^2;
end
srcnode = find(distance == min(distance));
%**************************************************************************

%assemble M,K
for ii = 1:1:num_elem
    % the nodes number
    i  = elements(ii,1);
    j  = elements(ii,2);
    k  = elements(ii,3);
    
    %the coordinate of the nodes
    xi = nodes(i,1);
    yi = nodes(i,2);
    
    xj = nodes(j,1);
    yj = nodes(j,2);
    
    xk = nodes(k,1);
    yk = nodes(k,2);
    
    %to calculate the grid distance
    grid_distance(ii,:) = pdist([xi,yi;xj,yj;xk,yk]);
    
    
    %the area of the triangle
    smatrix = [1,xi,yi;1,xj,yj;1,xk,yk];
    s       = 0.5*abs(det(smatrix));
       
    %M matrix
    M(i,i)  = s*(1/6)  + M(i,i);
    M(i,j)  = s*(1/12) + M(i,j);
    M(i,k)  = s*(1/12) + M(i,k);
    
    M(j,i)  = s*(1/12) + M(j,i);
    M(j,j)  = s*(1/6)  + M(j,j);
    M(j,k)  = s*(1/12) + M(j,k);
    
    M(k,i)  = s*(1/12) + M(k,i);
    M(k,j)  = s*(1/12) + M(k,j);
    M(k,k)  = s*(1/6)  + M(k,k);
    
    a = xi - xk;
    b = xj - xk;  
    c = yi - yk;
    d = yj - yk;
    cons = v^2/(4*s);

    %K matrix
    K(i,i) = K(i,i) + cons*(d*d+b*b);
    K(i,j) = K(i,j) + cons*(-d*c-b*a);
    K(i,k) = K(i,k) + cons*(d*(c-d)+b*(a-b));
    
    K(j,i) = K(j,i) + cons*(-c*d-a*b);
    K(j,j) = K(j,j) + cons*(c*c+a*a);
    K(j,k) = K(j,k) + cons*(c*(d-c)+a*(b-a));
    
    K(k,i) = K(k,i) + cons*(d*(c-d)+b*(a-b));
    K(k,j) = K(k,j) + cons*(c*(d-c)+a*(b-a));
    K(k,k) = K(k,k) + cons*((c-d)^2+(b-a)^2);
   
end
%output the maximum and minimum grid distance
maxgrid  =  max(max(grid_distance));
mingrid  =  min(min(grid_distance));
meangrid = 0.5*(maxgrid + mingrid);

%the time series of source wavelet.
% f0  = v/(ppw*h);%main frequency(Hz).
f0  = 25;
t   = 0:dt:(nt-1)*dt;
t0  = 1/f0;
amp = 1.0e7;
src = amp*(1 - 2.*(pi*f0.*(t - t0)).^2).*exp( - (pi*f0.*(t - t0)).^2);
%**************************************************************************
%calculate 
%sparse
temp1  = full((dt)^2*(II/M));

% %for gpu acceleration
% temp1g = gpuArray(temp1);wxw20230914modify
temp1g = temp1;

formatSpec = 'The calculating time node is: it = %d.\n';
for it =1:1:nt
    fprintf(formatSpec,it);
 
    F(srcnode,1)   = src(1,it); %load
    
    %sparse
    temp2 = full(F-K'*P);

    % gpu acceleration
    % temp2g = gpuArray(temp2); wxw20230914modify
    temp2g = temp2;
    temp = temp1g*temp2g;
    temp = gather(temp);
    Pnew  = temp + 2*P - Pold;
        
    Pold        = P;
    P           = Pnew;
    RES(:,it)   = Pnew;
 end
% post-processing
filename  = 'FEM2D.gif';
x        = nodes(:,1);
y        = nodes(:,2);
[xx, yy] = meshgrid(x_min:x_max,y_min:y_max);
zz = zeros(x_max+1, y_max+1, nt);
figure('color','w');
figure(1);
for ii = 1:1:nt
    time = ii*dt;%the present time(s)
    zz(:,:,ii) = griddata(x, y, RES(:,ii), xx, yy);
    imagesc(xx(1,:),yy(:,1),zz(:,:,ii));
    caxis([-0.3 0.3]);
    shading interp;
    colorbar;
    title(sprintf( 'Time Step = %d ; Time = %0.3f s ',ii,time));
    xlabel('X(m)');
    ylabel('Y(m)');
    set(gca,'FontName','Times New Roman','FontSize',15);
    drawnow; 
    FF = getframe(gcf);
    I = frame2im(FF); 
    [I, map] = rgb2ind(I, 256); 
    if ii == 1
        imwrite(I,map, filename,'gif', 'Loopcount',inf,'DelayTime',0.1);
    else
        imwrite(I,map, filename,'gif','WriteMode','append','DelayTime',0.1);
    end
end