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

 ``` 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) ```