[7e9223]: devel / example / Darcy-Stokes / P2P0 / Darcy_Stokes.m Maximize Restore History

Download this file

Darcy_Stokes.m    88 lines (60 with data), 2.3 kB

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
pkg load msh
pkg load fem-fenics
fem_init_env ();
fem_create_all ('DS');
fem_create_functional ('Err_v');
fem_create_functional ('Err_eps');
fem_create_functional ('Err_p');
uex = @(x, y) [-2*pi*sin(pi*x)*sin(pi*x)*cos(pi*y)*sin(pi*y);...
2*pi*sin(pi*x)*cos(pi*x)*sin(pi*y)*sin(pi*y)];
pex = @(x, y) sin(pi*x);
gradp = @(x, y) [pi*cos(pi*x); 0];
laplu = @(x, y) [-4*pi^3*(1 - 4*sin(pi*x)*sin(pi*x))*sin(pi*y)*cos(pi*y);...
4*pi^3*(1 - 4*sin(pi*y)*sin(pi*y))*sin(pi*x)*cos(pi*x)];
eps_vector = [1 2^-2 2^-4 2^-8 0]';
nel_vector = [4 8 16 32 64]';
h_vector = 1./nel_vector;
for ii = 1:numel(nel_vector)
nel = nel_vector(ii);
x = y = linspace (0, 1, nel + 1);
msho = msh2m_structured_mesh (x, y, 1, 1:4);
mshd = Mesh (msho);
v = Expression ('v', uex);
p = Expression ('p', pex);
V = DS_FunctionSpace (mshd);
V0 = SubSpace (V, 0);
bc = DirichletBC (V0, @(x,y) [0, 0], 1:4);
for jj = 1:numel(eps_vector)
epsi = eps_vector(jj);
ff = @(x, y) uex (x, y) - epsi * epsi * laplu (x, y) - gradp (x, y);
f = Expression ('f', ff);
g = Expression ('g', @(x,y) 0);
ep = Expression ('ep', @(x, y) epsi);
a = DS_BilinearForm (V, ep);
L = DS_LinearForm (V, f, g);
[A, b] = assemble_system (a, L, bc);
sol = A \ b;
func = Function ('uh', V, sol);
vh = Function ('vh', func, 0);
ph = Function ('ph', func, 1);
M = Err_v_Functional (mshd, v, vh);
err_v_L2(ii,jj) = sqrt (assemble (M));
const = fem_eval (ph, [0, 0]);
pex = @(x, y) sin(pi*x) + const;
p = Expression ('p', pex);
M = Err_p_Functional (mshd, p, ph);
err_p_L2(ii,jj) = sqrt (assemble (M));
M = Err_eps_Functional (mshd, v, vh, ep);
err_v_eps(ii,jj) = sqrt (assemble (M));
end
end
err_v_L2
err_p_L2
err_v_eps
pause ()
rate_v_L2 = (log (err_v_L2 (1:end-1,:)) - log (err_v_L2 (2:end,:)))./((log (h_vector (1:end-1)) - log (h_vector (2:end))));
conv_v_L2 = mean (rate_v_L2)
rate_p_L2 = (log (err_p_L2 (1:end-1,:)) - log (err_p_L2 (2:end,:)))./((log (h_vector (1:end-1)) - log (h_vector (2:end))));
conv_p_L2 = mean (rate_p_L2)
rate_v_eps = (log (err_v_eps (1:end-1,:)) - log (err_v_eps (2:end,:)))./((log (h_vector (1:end-1)) - log (h_vector (2:end))));
conv_v_eps = mean (rate_v_eps)