用节点DG求解二维热方程包括创建了一个特殊的集成电路用于测试。

clear all; close all; clc;

%% Preprocessing

% Parameters
 cfl = 0.2;
   N = 50;
tEnd = 0.05;

% build grid
x = linspace(-1,1,N+1); dx = 1/N;
y = linspace(-1,1,N+1); dy = 1/N;
[x,y] = meshgrid(x,y);

% set thermal diffusivity coef
k=0.1;

% Set IC
u0 = IC2d(x,y,2);

% set BCs
u0(1,1:N) = 0.1; 
u0(N,1:N) = 0.1;  
u0(1:N,1) = 0.1;  
u0(1:N,N) = 0.1;

%% Solver Loop

% load initial conditions 
u=u0; frame=0; iter=0; t=0; loop=1;

% set time step
dt=cfl*dx^2/k;

% initilize arrays
Laplace_u=zeros(size(u0));

while t < tEnd
    % time and counter
    t=t+dt; iter=iter+1;
    
    % load old step
    uo=u;
    
    % compute laplace operator
    for i = 2:N-1
        for j = 2:N-1
            Laplace_u(i,j) = (uo(i+1,j)-2*uo(i,j)+uo(i-1,j))/dx^2 + ...
                             (uo(i,j+1)-2*uo(i,j)+uo(i,j-1))/dy^2 ;
        end
    end
    
    % next time step
    u = uo + k*dt*Laplace_u;
    
    % get frame
    if (mod(iter,1)==0) % displays movie frame every 50 time steps
        frame=frame+1;
        surf(x,y,u); view(-45,45);
        h=gca;
        get(h,'FontSize');
        set(h,'FontSize',12);
        xlabel('X','fontSize',12);
        ylabel('Y','fontSize',12);
        colorbar('East');
        title('Heat Diffusion Equation','fontsize',12);
        fh = figure(1);
        set(fh, 'color', 'white');
        F=getframe;
    end
end

%% Postprocessing
% make a movie with the frames
movie(F,frame,1);

function u0 = IC2d(x,y,option)


switch option
    case 1 % 4 quadrants in Domain [-1,1]^2
        % Preallocate u0,
        u0 = zeros(size(x));

        % Initial Condition for our 2D domain
        u0(x >0 & y >0) = 0.50; % region 1
        u0(x<=0 & y >0) = 0.70; % region 2
        u0(x<=0 & y<=0) = 0.10; % region 3
        u0(x >0 & y<=0) = 0.90; % region 4
        
    case 2 % Square Jump
        % set center
        x0 = 0.0; y0 = 0.0;
        
        % parameters
        Dx = x(end)-x(1);
        Dy = y(end)-y(1);
        
        % Preallocate u0
        u0 = ones(size(x));
        
        % Parameters of region
        x1 = x0+Dx/10; x2 = x0-Dx/10;
        y1 = y0+Dy/10; y2 = y0-Dy/10;
        
        % Build Jump
        u0(x>x1) = 0.1;     u0(x<x2) = 0.1;
        u0(y>y1) = 0.1;     u0(y<y2) = 0.1;
        
    case 3 % Sine*Cosine 2-D in Domain [-1,1]^2
        u0 = sin(pi*x).*cos(pi*y);
        
    case 4 % Gaussian Jump
        % set center
        x0 = 0.0; y0 = 0.0; 
        
        % Gaussian
        u0 = 0.1 + 0.5*exp(-20*((x-x0).^2+(y-y0).^2));
        
    case 5 % Cylindrical Jump
        % set center
        x0 = 0.0; y0 = 0.0;
        
        % radious
        r = 0.5;
        
        % Gaussian
        u0 = 0.1*ones(size(x));
        u0(sqrt((x+x0).^2+(y+y0).^2)<r) = 1.0;
        
    case 6 % rectangle in x direction
        % set center
        x0 = 0.0;
        
        % parameters
        Dx = x(end)-x(1);
                
        % Preallocate u0
        u0 = ones(size(x));
        
        % Parameters of region
        x1 = x0+Dx/12; x2 = x0-Dx/12;
                
        % Build Jump
        u0(x>x1) = 0.1;     u0(x<x2) = 0.1;
                
  case 7 % rectangle in y direction
        % set center
        y0 = 0.0;
        
        % parameters
        Dy = y(end)-y(1);
        
        % Preallocate u0
        u0 = ones(size(x));
        
        % Parameters of region
        y1 = y0+Dy/12; y2 = y0-Dy/12;
        
        % Build Jump
        u0(y>y1) = 0.1;     u0(y<y2) = 0.1;
    case 8 % Riemann for range [-1,1]
        x1=x<0.0; x2=x>=0.0; 
        
        % Riemann IC
        u0=1*x1+0.1*x2;
        
    otherwise 
        error('case not listed :P')
end