Mercurial > repos > public > sbplib
comparison +sbp/+implementations/d2_variable_periodic_2.m @ 812:6b83dcb46f54 feature/grids
Merge with feature/poroelastic
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Fri, 27 Jul 2018 10:31:51 -0700 |
| parents | bbf303c1f0cf |
| children |
comparison
equal
deleted
inserted
replaced
| 798:e76321b89c1e | 812:6b83dcb46f54 |
|---|---|
| 1 function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_2(m,h) | |
| 2 % m = number of unique grid points, i.e. h = L/m; | |
| 3 | |
| 4 if(m<3) | |
| 5 error(['Operator requires at least ' num2str(3) ' grid points']); | |
| 6 end | |
| 7 | |
| 8 % Norm | |
| 9 Hv = ones(m,1); | |
| 10 Hv = h*Hv; | |
| 11 H = spdiag(Hv, 0); | |
| 12 HI = spdiag(1./Hv, 0); | |
| 13 | |
| 14 | |
| 15 % Dummy boundary operators | |
| 16 e_l = sparse(m,1); | |
| 17 e_r = rot90(e_l, 2); | |
| 18 | |
| 19 d1_l = sparse(m,1); | |
| 20 d1_r = -rot90(d1_l, 2); | |
| 21 | |
| 22 % D1 operator | |
| 23 diags = -1:1; | |
| 24 stencil = [-1/2 0 1/2]; | |
| 25 D1 = stripeMatrixPeriodic(stencil, diags, m); | |
| 26 D1 = D1/h; | |
| 27 | |
| 28 scheme_width = 3; | |
| 29 scheme_radius = (scheme_width-1)/2; | |
| 30 | |
| 31 r = 1:m; | |
| 32 offset = scheme_width; | |
| 33 r = r + offset; | |
| 34 | |
| 35 function D2 = D2_fun(c) | |
| 36 c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ]; | |
| 37 | |
| 38 Mm1 = -c(r-1)/2 - c(r)/2; | |
| 39 M0 = c(r-1)/2 + c(r) + c(r+1)/2; | |
| 40 Mp1 = -c(r)/2 - c(r+1)/2; | |
| 41 | |
| 42 vals = [Mm1,M0,Mp1]; | |
| 43 diags = -scheme_radius : scheme_radius; | |
| 44 M = spdiagsPeriodic(vals,diags); | |
| 45 | |
| 46 M=M/h; | |
| 47 D2=HI*(-M ); | |
| 48 end | |
| 49 D2 = @D2_fun; | |
| 50 end |