Mercurial > repos > public > sbplib

comparison +scheme/Beam.m @ 427:a613960a157b feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merged with feature/beams
author Ylva Rydin <ylva.rydin@telia.com>
date Thu, 26 Jan 2017 15:59:25 +0100
parents 25ae8333ce2b
children dfb27a7e801f
comparison
equal deleted inserted replaced
426:29944ea7674b 427:a613960a157b
1 classdef Beam < scheme.Scheme
2 properties
3 order % Order accuracy for the approximation
4 grid
5
6 D % non-stabalized scheme operator
7 alpha
8
9 h
10 H % Discrete norm
11 Hi
12
13 e_l, e_r
14 d1_l, d1_r
15 d2_l, d2_r
16 d3_l, d3_r
17 gamm
18 delt
19 alphaII
20 alphaIII
21
22 opt
23 end
24
25 methods
26 function obj = Beam(grid, order, alpha, opsGen, opt)
27 default_arg('alpha', -1);
28
29 % default_arg('opsGen', @sbp.D4);
30 default_arg('opsGen', @sbp.D4Variable); % Supposed to be better
31
32 opt_default.interface_l.tuning = 1.1;
33 opt_default.interface_l.tau = [];
34 opt_default.interface_l.sig = [];
35 opt_default.interface_r.tuning = 1.1;
36 opt_default.interface_r.tau = [];
37 opt_default.interface_r.sig = [];
38 default_struct('opt', opt_default);
39
40 if ~isa(grid, 'grid.Cartesian') || grid.D() ~= 1
41 error('Grid must be 1d cartesian');
42 end
43
44 obj.grid = grid;
45 obj.order = order;
46 obj.alpha = alpha;
47
48 m = grid.m;
49 h = grid.scaling();
50
51 x_lim = {grid.x{1}(1), grid.x{1}(end)};
52 ops = opsGen(m, x_lim, order);
53
54 D4 = ops.D4;
55 obj.H = ops.H;
56 obj.Hi = ops.HI;
57 obj.e_l = ops.e_l;
58 obj.e_r = ops.e_r;
59 obj.d1_l = ops.d1_l;
60 obj.d1_r = ops.d1_r;
61 obj.d2_l = ops.d2_l;
62 obj.d2_r = ops.d2_r;
63 obj.d3_l = ops.d3_l;
64 obj.d3_r = ops.d3_r;
65
66 obj.D = alpha*D4;
67
68 alphaII = ops.borrowing.N.S2/2;
69 alphaIII = ops.borrowing.N.S3/2;
70
71 obj.gamm = h*alphaII;
72 obj.delt = h^3*alphaIII;
73 obj.alphaII = alphaII;
74 obj.alphaIII = alphaIII;
75 obj.h = h;
76 obj.opt = opt;
77 end
78
79
80 % Closure functions return the opertors applied to the own doamin to close the boundary
81 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
82 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
83 % type is a string specifying the type of boundary condition if there are several.
84 % neighbour_scheme is an instance of Scheme that should be interfaced to.
85 % neighbour_boundary is a string specifying which boundary to interface to.
86 function [closure, penalty] = boundary_condition(obj,boundary,type)
87 default_arg('type','dn');
88
89 [e, d1, d2, d3, s] = obj.get_boundary_ops(boundary);
90 gamm = obj.gamm;
91 delt = obj.delt;
92
93 switch type
94 case {'dn', 'clamped'} % Dirichlet-neumann boundary condition
95 alpha = obj.alpha;
96
97 % tau1 < -alpha^2/gamma
98 % tuning = 2;
99 tuning = 1.1;
100
101 tau1 = tuning * alpha/delt;
102 tau4 = s*alpha;
103
104 sig2 = tuning * alpha/gamm;
105 sig3 = -s*alpha;
106
107 tau = tau1*e+tau4*d3;
108 sig = sig2*d1+sig3*d2;
109
110 closure = obj.Hi*(tau*e' + sig*d1');
111
112 penalty{1} = -obj.Hi*tau;
113 penalty{2} = -obj.Hi*sig;
114
115
116 case {'free'}
117 a = obj.alpha;
118
119 tau = s*a*d1;
120 sig = -s*a*e;
121
122 closure = obj.Hi*(tau*d2' + sig*d3');
123 penalty{1} = -obj.Hi*tau;
124 penalty{1} = -obj.Hi*sig;
125
126
127 otherwise % Unknown, boundary condition
128 error('No such boundary condition: type = %s',type);
129 end
130 end
131
132 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
133 % u denotes the solution in the own domain
134 % v denotes the solution in the neighbour domain
135 [e_u,d1_u,d2_u,d3_u,s_u] = obj.get_boundary_ops(boundary);
136 [e_v,d1_v,d2_v,d3_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
137
138
139 alpha_u = obj.alpha;
140 alpha_v = neighbour_scheme.alpha;
141
142
143 switch boundary
144 case 'l'
145 interface_opt = obj.opt.interface_l;
146 case 'r'
147 interface_opt = obj.opt.interface_r;
148 end
149
150
151 if isempty(interface_opt.tau) && isempty(interface_opt.sig)
152 gamm_u = obj.gamm;
153 delt_u = obj.delt;
154
155 gamm_v = neighbour_scheme.gamm;
156 delt_v = neighbour_scheme.delt;
157
158 tuning = interface_opt.tuning;
159
160 tau1 = ((alpha_u/2)/delt_u + (alpha_v/2)/delt_v)/2*tuning;
161 sig2 = ((alpha_u/2)/gamm_u + (alpha_v/2)/gamm_v)/2*tuning;
162 else
163 h_u = obj.h;
164 h_v = neighbour_scheme.h;
165
166 switch neighbour_boundary
167 case 'l'
168 neighbour_interface_opt = neighbour_scheme.opt.interface_l;
169 case 'r'
170 neighbour_interface_opt = neighbour_scheme.opt.interface_r;
171 end
172
173 tau_u = interface_opt.tau;
174 sig_u = interface_opt.sig;
175 tau_v = neighbour_interface_opt.tau;
176 sig_v = neighbour_interface_opt.sig;
177
178 tau1 = tau_u/h_u^3 + tau_v/h_v^3;
179 sig2 = sig_u/h_u + sig_v/h_v;
180 end
181
182 tau4 = s_u*alpha_u/2;
183 sig3 = -s_u*alpha_u/2;
184 phi2 = s_u*1/2;
185 psi1 = -s_u*1/2;
186
187 tau = tau1*e_u + tau4*d3_u;
188 sig = sig2*d1_u + sig3*d2_u ;
189 phi = phi2*d1_u ;
190 psi = psi1*e_u ;
191
192 closure = obj.Hi*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u');
193 penalty = -obj.Hi*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v');
194 end
195
196 % Returns the boundary ops and sign for the boundary specified by the string boundary.
197 % The right boundary is considered the positive boundary
198 function [e, d1, d2, d3, s] = get_boundary_ops(obj,boundary)
199 switch boundary
200 case 'l'
201 e = obj.e_l;
202 d1 = obj.d1_l;
203 d2 = obj.d2_l;
204 d3 = obj.d3_l;
205 s = -1;
206 case 'r'
207 e = obj.e_r;
208 d1 = obj.d1_r;
209 d2 = obj.d2_r;
210 d3 = obj.d3_r;
211 s = 1;
212 otherwise
213 error('No such boundary: boundary = %s',boundary);
214 end
215 end
216
217 function N = size(obj)
218 N = obj.grid.N;
219 end
220
221 end
222 end