Mercurial > repos > public > sbplib

comparison +rv/ResidualViscosity.m @ 1186:3364a51f0d9e feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add 1d shapiro filter function to ResidualViscosity
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 05 Jul 2019 19:08:37 +0200
parents beecb580c5bf
children 7173b6fd4063
comparison
equal deleted inserted replaced
1185:abb1b3ab8c23 1186:3364a51f0d9e
30 obj.Mres = obj.Cres*obj.h^2; 30 obj.Mres = obj.Cres*obj.h^2;
31 obj.Mfirst = obj.Cmax*obj.h; 31 obj.Mfirst = obj.Cmax*obj.h;
32 switch postProcess 32 switch postProcess
33 case 'none' 33 case 'none'
34 obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v); 34 obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v);
35 case 'filter'
36 F = obj.shapiroFilter(obj.grid, order);
37 obj.fRes = @(v,dvdt) F*obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v);
35 case 'maximum neighbors' 38 case 'maximum neighbors'
36 switch grid.D() 39 switch grid.D()
37 case 1 40 case 1
38 obj.fRes = @(v,dvdt) movmax(obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v),3); 41 obj.fRes = @(v,dvdt) movmax(obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v),3);
39 case 2 42 case 2
71 umax = max(movmax(uMatrix,3,1),movmax(uMatrix,3,2)); 74 umax = max(movmax(uMatrix,3,1),movmax(uMatrix,3,2));
72 umin = min(movmin(uMatrix,3,1),movmin(uMatrix,3,2)); 75 umin = min(movmin(uMatrix,3,1),movmin(uMatrix,3,2));
73 minmaxDiff = umax - umin; 76 minmaxDiff = umax - umin;
74 minmaxDiff = reshape(minmaxDiff',g.N(),1); 77 minmaxDiff = reshape(minmaxDiff',g.N(),1);
75 end 78 end
79 function F = shapiroFilter(grid, order)
80 switch order
81 case 2
82 F = spdiags(repmat([1 2 1],grid.m(1),1), -1:1, grid.m(1), grid.m(1));
83 case 4
84 F = spdiags(repmat([-1 4 10 4 -1],grid.m(1),1), -2:2, grid.m(1), grid.m(1));
85 case 6
86 F = spdiags(repmat([1 -6 15 44 15 -6 1],grid.m(1),1), -3:3, grid.m(1), grid.m(1));
87 case 8
88 F = spdiags(repmat([-1 8 -28 56 186 56 -28 8 -1],grid.m(1),1), -4:4, grid.m(1), grid.m(1));
89 end
90 F = 1/2^order * F;
91 % Fx = spdiags(repmat(1/(mid+2)*[1 mid 1],g.m(1),1), -1:1, g.m(1), g.m(1));
92 % Fx(1,:) = 0;
93 % Fx(1,1:2) = 1/(mid+1)*[mid 1];
94 % Fx(end,:) = 0;
95 % Fx(end,end-1:end) = 1/(mid+1)*[1 mid];
96 % Fy = spdiags(repmat(1/(mid+2)*[1 mid 1],g.m(2),1), -1:1, g.m(2), g.m(2));
97 % Fy(1,:) = 0;
98 % Fy(1,1:2) = 1/(mid+1)*[mid 1];
99 % Fy(end,:) = 0;
100 % Fy(end,end-1:end) = 1/(mid+1)*[1 mid];
101 % Ix = speye(g.m(1));
102 % Iy = speye(g.m(1));
103 % F = kron(Fx, Iy) + kron(Ix, Fy);
104 end
76 end 105 end
106
77 end 107 end