Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_variable_hollow_2.m @ 1303:49e3870335ef feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make the hollow scheme generation more efficient by introducing the D2VariableHollow opSet
author Martin Almquist <malmquist@stanford.edu>
date Sat, 11 Jul 2020 06:54:15 -0700
parents
children b5907140c069
comparison
equal deleted inserted replaced
1302:a0d615bde7f8 1303:49e3870335ef
1 function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_hollow_2(m,h)
2
3 BP = 1;
4 if(m<2*BP)
5 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
6 end
7
8 % Norm
9 Hv = ones(m,1);
10 Hv(1) = 1/2;
11 Hv(m:m) = 1/2;
12 Hv = h*Hv;
13 H = spdiag(Hv, 0);
14 HI = spdiag(1./Hv, 0);
15
16
17 % Boundary operators
18 e_l = sparse(m,1);
19 e_l(1) = 1;
20 e_r = rot90(e_l, 2);
21
22 d1_l = sparse(m,1);
23 d1_l(1:3) = 1/h*[-3/2 2 -1/2];
24 d1_r = -rot90(d1_l, 2);
25
26 % D1 operator
27 diags = -1:1;
28 stencil = [-1/2 0 1/2];
29 D1 = stripeMatrix(stencil, diags, m);
30
31 D1(1,1)=-1;D1(1,2)=1;D1(m,m-1)=-1;D1(m,m)=1;
32 D1(m,m-1)=-1;D1(m,m)=1;
33 D1=D1/h;
34 %Q=H*D1 + 1/2*(e_1*e_1') - 1/2*(e_m*e_m');
35
36
37 M=sparse(m,m);
38 nBP = 2;
39
40 scheme_width = 3;
41 scheme_radius = (scheme_width-1)/2;
42
43 function D2 = D2_fun(c)
44
45 for r = 1+scheme_radius:nBP
46 Mm1 = -c(r-1)/2 - c(r)/2;
47 M0 = c(r-1)/2 + c(r) + c(r+1)/2;
48 Mp1 = -c(r)/2 - c(r+1)/2;
49
50 M(r, r-1) = Mm1;
51 M(r, r) = M0;
52 M(r, r+1) = Mp1;
53 end
54
55 for r = m-nBP+1:m-scheme_radius
56 Mm1 = -c(r-1)/2 - c(r)/2;
57 M0 = c(r-1)/2 + c(r) + c(r+1)/2;
58 Mp1 = -c(r)/2 - c(r+1)/2;
59
60 M(r, r-1) = Mm1;
61 M(r, r) = M0;
62 M(r, r+1) = Mp1;
63 end
64
65 M(1:2,1:2)=[c(1)/2 + c(2)/2 -c(1)/2 - c(2)/2; -c(1)/2 - c(2)/2 c(1)/2 + c(2) + c(3)/2;];
66 M(m-1:m,m-1:m)=[c(m-2)/2 + c(m-1) + c(m)/2 -c(m-1)/2 - c(m)/2; -c(m-1)/2 - c(m)/2 c(m-1)/2 + c(m)/2;];
67 M=M/h;
68
69 D2=HI*(-M-c(1)*e_l*d1_l'+c(m)*e_r*d1_r');
70 end
71 D2 = @D2_fun;
72 end