Mercurial > repos > public > sbplib
comparison +sbp/+implementations/d2_variable_periodic_2.m @ 681:7368affc8f78 feature/poroelastic
Add D2 variable periodic for second order.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Wed, 07 Feb 2018 15:42:50 -0800 |
| parents | |
| children | 5ccf6aaf6d6b |
comparison
equal
deleted
inserted
replaced
| 680:cd1a76c38565 | 681:7368affc8f78 |
|---|---|
| 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 BP = 1; | |
| 5 if(m<2*BP) | |
| 6 error(['Operator requires at least ' num2str(2*BP) ' grid points']); | |
| 7 end | |
| 8 | |
| 9 % Norm | |
| 10 Hv = ones(m,1); | |
| 11 Hv = h*Hv; | |
| 12 H = spdiag(Hv, 0); | |
| 13 HI = spdiag(1./Hv, 0); | |
| 14 | |
| 15 | |
| 16 % Dummy boundary operators | |
| 17 e_l = sparse(m,1); | |
| 18 e_r = rot90(e_l, 2); | |
| 19 | |
| 20 d1_l = sparse(m,1); | |
| 21 d1_r = -rot90(d1_l, 2); | |
| 22 | |
| 23 % D1 operator | |
| 24 diags = -1:1; | |
| 25 stencil = [-1/2 0 1/2]; | |
| 26 D1 = stripeMatrixPeriodic(stencil, diags, m); | |
| 27 D1 = D1/h; | |
| 28 | |
| 29 scheme_width = 3; | |
| 30 scheme_radius = (scheme_width-1)/2; | |
| 31 | |
| 32 r = 1:m; | |
| 33 offset = scheme_width; | |
| 34 r = r + offset; | |
| 35 | |
| 36 function D2 = D2_fun(c) | |
| 37 c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ]; | |
| 38 | |
| 39 Mm1 = -c(r-1)/2 - c(r)/2; | |
| 40 M0 = c(r-1)/2 + c(r) + c(r+1)/2; | |
| 41 Mp1 = -c(r)/2 - c(r+1)/2; | |
| 42 | |
| 43 vals = [Mm1,M0,Mp1]; | |
| 44 diags = -scheme_radius : scheme_radius; | |
| 45 M = spdiagsVariablePeriodic(vals,diags); | |
| 46 | |
| 47 M=M/h; | |
| 48 D2=HI*(-M ); | |
| 49 end | |
| 50 D2 = @D2_fun; | |
| 51 end |