43 lines
1.1 KiB
Matlab
43 lines
1.1 KiB
Matlab
clear; clc; close all;
|
|
p = 70;
|
|
%p=-70;
|
|
lambda = @(x) 1; %constant diffusion
|
|
f = @(x) 0; %homogeneous RHS
|
|
|
|
N = 30; % Initial mesh (uniform)
|
|
x = linspace(0,1,N)';
|
|
|
|
nIter = 8;
|
|
|
|
for it = 1:nIter
|
|
[K,F] = assemble_1D(x,lambda,f); %assemble
|
|
for i = 1:length(x)-1 %convection term added manually
|
|
h = x(i+1) - x(i);
|
|
Ke_conv = p/2 * [-1 1; -1 1];
|
|
K(i:i+1,i:i+1) = K(i:i+1,i:i+1) + Ke_conv;
|
|
end
|
|
K(1,:) = 0; K(1,1) = 1; F(1) = 0; % Dirichlet BC at x=0
|
|
K(end,:) = 0; K(end,end) = 1; F(end) = 1; % Dirichlet BC at x=1
|
|
|
|
u = K\F;
|
|
|
|
if it < nIter
|
|
Nn = length(x); %initialize indicator
|
|
eta = zeros(Nn,1);
|
|
for j = 2:Nn-1 %interior nodes only
|
|
ul = (u(j)-u(j-1))/(x(j)-x(j-1)); %left derivative
|
|
ur = (u(j+1)-u(j))/(x(j+1)-x(j)); %right derivative
|
|
%Jump in convection-diffusion flux
|
|
Jl = -ul + p*u(j);
|
|
Jr = -ur + p*u(j);
|
|
eta(j) = abs(Jr - Jl);
|
|
end
|
|
x = r_adapt(x, eta);
|
|
end
|
|
end
|
|
|
|
figure
|
|
plot(x,u,'-o','LineWidth',1.5)
|
|
xlabel('x'), ylabel('u')
|
|
title(['Ex6C, r-adaptivity, Péclet p = ',num2str(p)])
|
|
grid on
|