Mercurial > repos > public > sbplib

comparison +scheme/Beam.m @ 240:46256fffa329 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added some stuff for default structs. Improved interface for beam.
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 03 Aug 2016 12:30:27 +0200
parents ba56e0d621f2
children f4306afd6849
comparison
equal deleted inserted replaced
239:ba56e0d621f2 240:46256fffa329
13 d1_l, d1_r 13 d1_l, d1_r
14 d2_l, d2_r 14 d2_l, d2_r
15 d3_l, d3_r 15 d3_l, d3_r
16 gamm 16 gamm
17 delt 17 delt
18 interface_tuning 18
19 opt
19 end 20 end
20 21
21 methods 22 methods
22 function obj = Beam(grid, order, alpha, opsGen, interface_tuning, alphaII, alphaIII) 23 function obj = Beam(grid, order, alpha, opsGen, opt)
23 default_arg('alpha', -1); 24 default_arg('alpha', -1);
24 default_arg('interface_tuning', 1.1); 25
25 default_arg('alphaII', []) 26
26 default_arg('alphaIII', []) 27 opt_default.interface_l.tuning = 1.1;
28 opt_default.interface_l.tau = [];
29 opt_default.interface_l.sig = [];
30 opt_default.interface_r.tuning = 1.1;
31 opt_default.interface_r.tau = [];
32 opt_default.interface_r.sig = [];
33 default_struct('opt', opt_default);
34
27 35
28 % default_arg('opsGen', @sbp.Higher); 36 % default_arg('opsGen', @sbp.Higher);
29 default_arg('opsGen', @sbp.HigherCompatibleVariable); % Supposed to be better 37 default_arg('opsGen', @sbp.HigherCompatibleVariable); % Supposed to be better
30 38
31 if ~isa(grid, 'grid.Cartesian') || grid.D() ~= 1 39 if ~isa(grid, 'grid.Cartesian') || grid.D() ~= 1
55 obj.d3_l = sparse(ops.boundary.S3_1); 63 obj.d3_l = sparse(ops.boundary.S3_1);
56 obj.d3_r = sparse(ops.boundary.S3_m); 64 obj.d3_r = sparse(ops.boundary.S3_m);
57 65
58 obj.D = alpha*D4; 66 obj.D = alpha*D4;
59 67
60 if isempty(alphaII) && isempty(alphaIII) 68 alphaII = ops.borrowing.N.S2/2;
61 alphaII = ops.borrowing.N.S2/2; 69 alphaIII = ops.borrowing.N.S3/2;
62 alphaIII = ops.borrowing.N.S3/2;
63 end
64 70
65 obj.gamm = h*alphaII; 71 obj.gamm = h*alphaII;
66 obj.delt = h^3*alphaIII; 72 obj.delt = h^3*alphaIII;
67 obj.interface_tuning = interface_tuning; 73 obj.opt = opt;
68 end 74 end
69 75
70 76
71 % Closure functions return the opertors applied to the own doamin to close the boundary 77 % Closure functions return the opertors applied to the own doamin to close the boundary
72 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 78 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
111 % u denotes the solution in the own domain 117 % u denotes the solution in the own domain
112 % v denotes the solution in the neighbour domain 118 % v denotes the solution in the neighbour domain
113 [e_u,d1_u,d2_u,d3_u,s_u] = obj.get_boundary_ops(boundary); 119 [e_u,d1_u,d2_u,d3_u,s_u] = obj.get_boundary_ops(boundary);
114 [e_v,d1_v,d2_v,d3_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 120 [e_v,d1_v,d2_v,d3_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
115 121
116 gamm_u = obj.gamm;
117 delt_u = obj.delt;
118
119 gamm_v = neighbour_scheme.gamm;
120 delt_v = neighbour_scheme.delt;
121
122 tuning = obj.interface_tuning;
123
124 122
125 alpha_u = obj.alpha; 123 alpha_u = obj.alpha;
126 alpha_v = neighbour_scheme.alpha; 124 alpha_v = neighbour_scheme.alpha;
127 125
128 tau1 = ((alpha_u/2)/delt_u + (alpha_v/2)/delt_v)/2*tuning; 126 switch boundary
129 % tau1 = (alpha_u/2 + alpha_v/2)/(2*delt_u)*tuning; 127 case 'l'
128 interface_opt = obj.opt.interface_l;
129 case 'r'
130 interface_opt = obj.opt.interface_r;
131 end
132
133
134 if isempty(interface_opt.tau) && isempty(interface_opt.sig)
135 gamm_u = obj.gamm;
136 delt_u = obj.delt;
137
138 gamm_v = neighbour_scheme.gamm;
139 delt_v = neighbour_scheme.delt;
140
141 tuning = interface_opt.tuning;
142
143 tau1 = ((alpha_u/2)/delt_u + (alpha_v/2)/delt_v)/2*tuning;
144 sig2 = ((alpha_u/2)/gamm_u + (alpha_v/2)/gamm_v)/2*tuning;
145 else
146 tau1 = interface_opt.tau;
147 sig2 = interface_opt.sig;
148 end
149
130 tau4 = s_u*alpha_u/2; 150 tau4 = s_u*alpha_u/2;
131
132 sig2 = ((alpha_u/2)/gamm_u + (alpha_v/2)/gamm_v)/2*tuning;
133 sig3 = -s_u*alpha_u/2; 151 sig3 = -s_u*alpha_u/2;
134
135 phi2 = s_u*1/2; 152 phi2 = s_u*1/2;
136
137 psi1 = -s_u*1/2; 153 psi1 = -s_u*1/2;
138 154
139 tau = tau1*e_u + tau4*d3_u; 155 tau = tau1*e_u + tau4*d3_u;
140 sig = sig2*d1_u + sig3*d2_u ; 156 sig = sig2*d1_u + sig3*d2_u ;
141 phi = phi2*d1_u ; 157 phi = phi2*d1_u ;