From 644766399f696ddcec5066a2d69f2604d78a4004 Mon Sep 17 00:00:00 2001 From: Wxwei <12221108@zju.edu.cn> Date: Mon, 18 Sep 2023 10:24:55 +0000 Subject: [PATCH] =?UTF-8?q?=E4=B8=8A=E4=BC=A0=E6=96=87=E4=BB=B6=E8=87=B3?= =?UTF-8?q?=20'FEM=5Ftest/WaveEq2D'?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 有限元方法案例，求解二维声波传播过程 --- FEM_test/WaveEq2D/D2waveEq.m | 170 +++++++++++++++++++++++++++++++++++ 1 file changed, 170 insertions(+) create mode 100644 FEM_test/WaveEq2D/D2waveEq.m diff --git a/FEM_test/WaveEq2D/D2waveEq.m b/FEM_test/WaveEq2D/D2waveEq.m new file mode 100644 index 0000000..90a94da --- /dev/null +++ b/FEM_test/WaveEq2D/D2waveEq.m @@ -0,0 +1,170 @@ +%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 \ No newline at end of file