Mercurial > repos > public > sbplib

comparison +time/SBPInTime.m @ 422:38173ea263ed feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in default
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 07 Feb 2017 15:59:22 +0100
parents 9ff24a14f9ef e1db62d14835
children 8894e9c49e40
comparison
equal deleted inserted replaced
402:bc78157c89cb 422:38173ea263ed
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;
96 end 105 end
97 106
98 107
99 methods(Static) 108 methods(Static)
100 function N = smallestBlockSize(order,TYPE) 109 function N = smallestBlockSize(order,TYPE)
101 default_arg('TYPE','equidistant') 110 default_arg('TYPE','gauss')
102 111
103 switch TYPE 112 switch TYPE
104 113
105 case 'equidistant' 114 case 'equidistant'
106 switch order 115 switch order
151 case 12 160 case 12
152 N = 20; 161 N = 20;
153 otherwise 162 otherwise
154 error('Operator does not exist'); 163 error('Operator does not exist');
155 end 164 end
165 case 'gauss'
166 N = 4;
156 end 167 end
157 end 168 end
158 end 169 end
159 end 170 end