Mercurial > repos > public > sbplib

comparison +scheme/Beam.m @ 998:2b1b944deae1 feature/getBoundaryOp

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add getBoundaryOperator to all 1d schemes. Did not add getBoundaryQuadrature because it doesnt make sense in 1d (?)
author Martin Almquist <malmquist@stanford.edu>
date Sat, 12 Jan 2019 13:35:19 -0800
parents a35ed1d124d3
children 8d73fcdb07a5
comparison
equal deleted inserted replaced
997:78db023a7fe3 998:2b1b944deae1
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, d1, d2, d3] = obj.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, boundary);
90 s = obj.getBoundarySign(boundary);
90 gamm = obj.gamm; 91 gamm = obj.gamm;
91 delt = obj.delt; 92 delt = obj.delt;
92 93
93 94
94 % TODO: Can this be simplifed? Can I handle conditions on u on its own, u_x on its own ... 95 % TODO: Can this be simplifed? Can I handle conditions on u on its own, u_x on its own ...
171 end 172 end
172 173
173 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary, type) 174 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary, type)
174 % u denotes the solution in the own domain 175 % u denotes the solution in the own domain
175 % v denotes the solution in the neighbour domain 176 % v denotes the solution in the neighbour domain
176 [e_u,d1_u,d2_u,d3_u,s_u] = obj.get_boundary_ops(boundary); 177 [e_u, d1_u, d2_u, d3_u] = obj.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, boundary);
177 [e_v,d1_v,d2_v,d3_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 178 s_u = obj.getBoundarySign(boundary);
178 179
180 [e_v, d1_v, d2_v, d3_v] = neighbour_scheme.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, neighbour_boundary);
181 s_v = neighbour_scheme.getBoundarySign(neighbour_boundary);
179 182
180 alpha_u = obj.alpha; 183 alpha_u = obj.alpha;
181 alpha_v = neighbour_scheme.alpha; 184 alpha_v = neighbour_scheme.alpha;
182
183 185
184 switch boundary 186 switch boundary
185 case 'l' 187 case 'l'
186 interface_opt = obj.opt.interface_l; 188 interface_opt = obj.opt.interface_l;
187 case 'r' 189 case 'r'
232 234
233 closure = obj.Hi*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u'); 235 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'); 236 penalty = -obj.Hi*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v');
235 end 237 end
236 238
237 % Returns the boundary ops and sign for the boundary specified by the string boundary. 239 % Returns the boundary operator op for the boundary specified by the string boundary.
238 % The right boundary is considered the positive boundary 240 % op -- string or a cell array of strings
239 function [e, d1, d2, d3, s] = get_boundary_ops(obj,boundary) 241 % boundary -- string
242 function varargout = getBoundaryOperator(obj, op, boundary)
243
244 if ~ismember(boundary, {'l', 'r'})
245 error('No such boundary: boundary = %s',boundary);
246 end
247
248 if ~iscell(op)
249 op = {op};
250 end
251
252 for i = 1:numel(op)
253 switch op{i}
254 case 'e'
255 switch boundary
256 case 'l'
257 e = obj.e_l;
258 case 'r'
259 e = obj.e_r;
260 end
261 varargout{i} = e;
262
263 case 'd1'
264 switch boundary
265 case 'l'
266 d1 = obj.d1_l;
267 case 'r'
268 d1 = obj.d1_r;
269 end
270 varargout{i} = d1;
271 end
272
273 case 'd2'
274 switch boundary
275 case 'l'
276 d2 = obj.d2_l;
277 case 'r'
278 d2 = obj.d2_r;
279 end
280 varargout{i} = d2;
281 end
282
283 case 'd3'
284 switch boundary
285 case 'l'
286 d3 = obj.d3_l;
287 case 'r'
288 d3 = obj.d3_r;
289 end
290 varargout{i} = d3;
291 end
292 end
293 end
294
295 % Returns the boundary sign. The right boundary is considered the positive boundary
296 % boundary -- string
297 function s = getBoundarySign(obj, boundary)
240 switch boundary 298 switch boundary
241 case 'l' 299 case {'r'}
242 e = obj.e_l; 300 s = 1;
243 d1 = obj.d1_l; 301 case {'l'}
244 d2 = obj.d2_l;
245 d3 = obj.d3_l;
246 s = -1; 302 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 303 otherwise
254 error('No such boundary: boundary = %s',boundary); 304 error('No such boundary: boundary = %s',boundary);
255 end 305 end
256 end 306 end
257 307