Mercurial > repos > public > sbplib

comparison +scheme/Beam.m @ 1197:433c89bf19e0 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 07 Aug 2019 15:23:42 +0200
parents 0c504a21432d
children
comparison
equal deleted inserted replaced
1196:f6c571d8f22f 1197:433c89bf19e0
84 % neighbour_scheme is an instance of Scheme that should be interfaced to. 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. 85 % neighbour_boundary is a string specifying which boundary to interface to.
86 function [closure, penalty] = boundary_condition(obj,boundary,type) 86 function [closure, penalty] = boundary_condition(obj,boundary,type)
87 default_arg('type','dn'); 87 default_arg('type','dn');
88 88
89 [e, d1, d2, d3, s] = obj.get_boundary_ops(boundary); 89 e = obj.getBoundaryOperator('e', boundary);
90 d1 = obj.getBoundaryOperator('d1', boundary);
91 d2 = obj.getBoundaryOperator('d2', boundary);
92 d3 = obj.getBoundaryOperator('d3', boundary);
93 s = obj.getBoundarySign(boundary);
90 gamm = obj.gamm; 94 gamm = obj.gamm;
91 delt = obj.delt; 95 delt = obj.delt;
92 96
93 97
94 % TODO: Can this be simplifed? Can I handle conditions on u on its own, u_x on its own ... 98 % TODO: Can this be simplifed? Can I handle conditions on u on its own, u_x on its own ...
122 tau = s*a*d1; 126 tau = s*a*d1;
123 sig = -s*a*e; 127 sig = -s*a*e;
124 128
125 closure = obj.Hi*(tau*d2' + sig*d3'); 129 closure = obj.Hi*(tau*d2' + sig*d3');
126 penalty{1} = -obj.Hi*tau; 130 penalty{1} = -obj.Hi*tau;
127 penalty{1} = -obj.Hi*sig; 131 penalty{2} = -obj.Hi*sig;
128 132
129 case 'e' 133 case 'e'
130 alpha = obj.alpha; 134 alpha = obj.alpha;
131 tuning = 1.1; 135 tuning = 1.1;
132 136
171 end 175 end
172 176
173 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary, type) 177 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary, type)
174 % u denotes the solution in the own domain 178 % u denotes the solution in the own domain
175 % v denotes the solution in the neighbour domain 179 % v denotes the solution in the neighbour domain
176 [e_u,d1_u,d2_u,d3_u,s_u] = obj.get_boundary_ops(boundary); 180 e_u = obj.getBoundaryOperator('e', boundary);
177 [e_v,d1_v,d2_v,d3_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 181 d1_u = obj.getBoundaryOperator('d1', boundary);
178 182 d2_u = obj.getBoundaryOperator('d2', boundary);
183 d3_u = obj.getBoundaryOperator('d3', boundary);
184 s_u = obj.getBoundarySign(boundary);
185
186 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
187 d1_v = neighbour_scheme.getBoundaryOperator('d1', neighbour_boundary);
188 d2_v = neighbour_scheme.getBoundaryOperator('d2', neighbour_boundary);
189 d3_v = neighbour_scheme.getBoundaryOperator('d3', neighbour_boundary);
190 s_v = neighbour_scheme.getBoundarySign(neighbour_boundary);
179 191
180 alpha_u = obj.alpha; 192 alpha_u = obj.alpha;
181 alpha_v = neighbour_scheme.alpha; 193 alpha_v = neighbour_scheme.alpha;
182
183 194
184 switch boundary 195 switch boundary
185 case 'l' 196 case 'l'
186 interface_opt = obj.opt.interface_l; 197 interface_opt = obj.opt.interface_l;
187 case 'r' 198 case 'r'
232 243
233 closure = obj.Hi*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u'); 244 closure = obj.Hi*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u');
234 penalty = -obj.Hi*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v'); 245 penalty = -obj.Hi*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v');
235 end 246 end
236 247
237 % Returns the boundary ops and sign for the boundary specified by the string boundary. 248 % Returns the boundary operator op for the boundary specified by the string boundary.
238 % The right boundary is considered the positive boundary 249 % op -- string
239 function [e, d1, d2, d3, s] = get_boundary_ops(obj,boundary) 250 % boundary -- string
251 function o = getBoundaryOperator(obj, op, boundary)
252 assertIsMember(op, {'e', 'd1', 'd2', 'd3'})
253 assertIsMember(boundary, {'l', 'r'})
254
255 o = obj.([op, '_', boundary]);
256 end
257
258 % Returns square boundary quadrature matrix, of dimension
259 % corresponding to the number of boundary points
260 %
261 % boundary -- string
262 % Note: for 1d diffOps, the boundary quadrature is the scalar 1.
263 function H_b = getBoundaryQuadrature(obj, boundary)
264 assertIsMember(boundary, {'l', 'r'})
265
266 H_b = 1;
267 end
268
269 % Returns the boundary sign. The right boundary is considered the positive boundary
270 % boundary -- string
271 function s = getBoundarySign(obj, boundary)
272 assertIsMember(boundary, {'l', 'r'})
273
240 switch boundary 274 switch boundary
241 case 'l' 275 case {'r'}
242 e = obj.e_l; 276 s = 1;
243 d1 = obj.d1_l; 277 case {'l'}
244 d2 = obj.d2_l;
245 d3 = obj.d3_l;
246 s = -1; 278 s = -1;
247 case 'r'
248 e = obj.e_r;
249 d1 = obj.d1_r;
250 d2 = obj.d2_r;
251 d3 = obj.d3_r;
252 s = 1;
253 otherwise
254 error('No such boundary: boundary = %s',boundary);
255 end 279 end
256 end 280 end
257 281
258 function N = size(obj) 282 function N = size(obj)
259 N = obj.grid.N; 283 N = obj.grid.N;