Mercurial > repos > public > sbplib
comparison +sbp/D2VariableCompatibleHollow.m @ 1303:49e3870335ef feature/poroelastic
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 classdef D2VariableCompatibleHollow < sbp.OpSet | |
| 2 properties | |
| 3 D1 % SBP operator approximating first derivative | |
| 4 H % Norm matrix | |
| 5 HI % H^-1 | |
| 6 Q % Skew-symmetric matrix | |
| 7 e_l % Left boundary operator | |
| 8 e_r % Right boundary operator | |
| 9 D2 % SBP operator for second derivative | |
| 10 M % Norm matrix, second derivative | |
| 11 d1_l % Left boundary first derivative | |
| 12 d1_r % Right boundary first derivative | |
| 13 m % Number of grid points. | |
| 14 h % Step size | |
| 15 x % grid | |
| 16 borrowing % Struct with borrowing limits for different norm matrices | |
| 17 end | |
| 18 | |
| 19 methods | |
| 20 function obj = D2VariableCompatibleHollow(m,lim,order) | |
| 21 | |
| 22 x_l = lim{1}; | |
| 23 x_r = lim{2}; | |
| 24 L = x_r-x_l; | |
| 25 obj.h = L/(m-1); | |
| 26 obj.x = linspace(x_l,x_r,m)'; | |
| 27 | |
| 28 switch order | |
| 29 | |
| 30 case 6 | |
| 31 | |
| 32 [obj.H, obj.HI, obj.D1, D2, ... | |
| 33 ~, obj.e_l, obj.e_r, ~, ~, ~, ~, ~,... | |
| 34 d1_l, d1_r] = ... | |
| 35 sbp.implementations.d4_variable_hollow_6(m, obj.h); | |
| 36 | |
| 37 case 4 | |
| 38 [obj.H, obj.HI, obj.D1, D2, obj.e_l,... | |
| 39 obj.e_r, d1_l, d1_r] = ... | |
| 40 sbp.implementations.d2_variable_hollow_4(m,obj.h); | |
| 41 case 2 | |
| 42 [obj.H, obj.HI, obj.D1, D2, obj.e_l,... | |
| 43 obj.e_r, d1_l, d1_r] = ... | |
| 44 sbp.implementations.d2_variable_hollow_2(m,obj.h); | |
| 45 | |
| 46 otherwise | |
| 47 error('Invalid operator order %d.',order); | |
| 48 end | |
| 49 obj.borrowing.H11 = obj.H(1,1)/obj.h; % First element in H/h, | |
| 50 obj.borrowing.M.d1 = obj.H(1,1)/obj.h; % First element in H/h is borrowing also for M | |
| 51 obj.borrowing.R.delta_D = inf; | |
| 52 obj.m = m; | |
| 53 obj.M = []; | |
| 54 | |
| 55 | |
| 56 D1 = obj.D1; | |
| 57 e_r = obj.e_r; | |
| 58 e_l = obj.e_l; | |
| 59 | |
| 60 % D2 = Hinv * (-M + br*er*d1r^T - bl*el*d1l^T); | |
| 61 % Replace d1' by e'*D1 in D2. | |
| 62 D2_compatible = @(b) D2(b) - obj.HI*(b(m)*e_r*d1_r' - b(m)*e_r*e_r'*D1) ... | |
| 63 + obj.HI*(b(1)*e_l*d1_l' - b(1)*e_l*e_l'*D1); | |
| 64 | |
| 65 obj.D2 = D2_compatible; | |
| 66 obj.d1_l = (e_l'*D1)'; | |
| 67 obj.d1_r = (e_r'*D1)'; | |
| 68 | |
| 69 end | |
| 70 function str = string(obj) | |
| 71 str = [class(obj) '_' num2str(obj.order)]; | |
| 72 end | |
| 73 end | |
| 74 | |
| 75 | |
| 76 end | |
| 77 | |
| 78 | |
| 79 | |
| 80 | |
| 81 |