Mercurial > repos > public > sbplib
comparison +sbp/+implementations/d2_variable_hollow_2.m @ 1331:60c875c18de3 feature/D2_boundary_opt
Merge with feature/poroelastic for Elastic schemes
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Thu, 10 Mar 2022 16:54:26 +0100 |
| parents | b5907140c069 |
| children |
comparison
equal
deleted
inserted
replaced
| 1330:855871e0b852 | 1331:60c875c18de3 |
|---|---|
| 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 nBP = 2; | |
| 38 M = sparse(m,m); | |
| 39 coeffs = load('sbplib/+sbp/+implementations/coeffs_d2_variable_2.mat'); | |
| 40 | |
| 41 function D2 = D2_fun(c) | |
| 42 M_l = zeros(nBP, coeffs.nBPC); | |
| 43 M_r = zeros(nBP, coeffs.nBPC); | |
| 44 | |
| 45 for i=1:coeffs.nBPC | |
| 46 M_l = M_l + coeffs.C_l{i}*c(i); | |
| 47 M_r = M_r + coeffs.C_r{i}*c(m-coeffs.nBPC+i); | |
| 48 end | |
| 49 | |
| 50 M(1:nBP, 1:coeffs.nBPC) = M_l; | |
| 51 M(m-nBP+1:m, m-coeffs.nBPC+1:m) = M_r; | |
| 52 | |
| 53 D2 = M/h^2; | |
| 54 end | |
| 55 D2 = @D2_fun; | |
| 56 end |