Mercurial > repos > public > sbplib

comparison +time/SBPInTime.m @ 411:16a76db8c279

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged in feature/SBPInTimeGauss (pull request #6) Feature/sbpintimegauss
author Jonatan Werpers <jonatan.werpers@it.uu.se>
date Tue, 07 Feb 2017 12:54:52 +0000
parents e1db62d14835
children 38173ea263ed b5e5b195da1e
comparison
equal deleted inserted replaced
404:d6d27fdc342a 411:16a76db8c279
22 nodes 22 nodes
23 end 23 end
24 24
25 methods 25 methods
26 function obj = SBPInTime(A, f, k, t0, v0, TYPE, order, blockSize) 26 function obj = SBPInTime(A, f, k, t0, v0, TYPE, order, blockSize)
27 default_arg('TYPE','minimal'); 27
28 default_arg('order', 8); 28 default_arg('TYPE','gauss');
29 default_arg('blockSize',time.SBPInTime.smallestBlockSize(order,TYPE)); 29
30 if(strcmp(TYPE,'gauss'))
31 default_arg('order',4)
32 default_arg('blockSize',4)
33 else
34 default_arg('order', 8);
35 default_arg('blockSize',time.SBPInTime.smallestBlockSize(order,TYPE));
36 end
30 37
31 obj.A = A; 38 obj.A = A;
32 obj.f = f; 39 obj.f = f;
33 obj.k_local = k/(blockSize-1); 40 obj.k_local = k/(blockSize-1);
34 obj.k = k; 41 obj.k = k;
43 ops = sbp.D2Standard(blockSize,{0,obj.k},order); 50 ops = sbp.D2Standard(blockSize,{0,obj.k},order);
44 case 'optimal' 51 case 'optimal'
45 ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order); 52 ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order);
46 case 'minimal' 53 case 'minimal'
47 ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal'); 54 ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal');
55 case 'gauss'
56 ops = sbp.D1Gauss(blockSize,{0,obj.k});
48 end 57 end
49 58
50 D1 = ops.D1; 59 D1 = ops.D1;
51 HI = ops.HI; 60 HI = ops.HI;
52 e_l = ops.e_l; 61 e_l = ops.e_l;
74 % LU factorization 83 % LU factorization
75 [obj.L,obj.U,obj.P,obj.Q] = lu(obj.M); 84 [obj.L,obj.U,obj.P,obj.Q] = lu(obj.M);
76 85
77 % Pretend that the initial condition is the last level 86 % Pretend that the initial condition is the last level
78 % of a previous step. 87 % of a previous step.
79 obj.v = obj.Et_r * v0; 88 obj.v = 1/(e_r'*e_r) * obj.Et_r * v0;
80 89
81 end 90 end
82 91
83 function [v,t] = getV(obj) 92 function [v,t] = getV(obj)
84 v = obj.Et_r' * obj.v; 93 v = obj.Et_r' * obj.v;
95 end 104 end
96 end 105 end
97 106
98 methods(Static) 107 methods(Static)
99 function N = smallestBlockSize(order,TYPE) 108 function N = smallestBlockSize(order,TYPE)
100 default_arg('TYPE','equidistant') 109 default_arg('TYPE','gauss')
101 110
102 switch TYPE 111 switch TYPE
103 112
104 case 'equidistant' 113 case 'equidistant'
105 switch order 114 switch order
150 case 12 159 case 12
151 N = 20; 160 N = 20;
152 otherwise 161 otherwise
153 error('Operator does not exist'); 162 error('Operator does not exist');
154 end 163 end
164 case 'gauss'
165 N = 4;
155 end 166 end
156 end 167 end
157 end 168 end
158 end 169 end