Mercurial > repos > public > sbplib

comparison +scheme/Beam.m @ 1108:5ec23b9bf360 feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Martin Almquist <malmquist@stanford.edu>
date Wed, 10 Apr 2019 11:00:27 -0700
parents 0c504a21432d
children
comparison
equal deleted inserted replaced
1087:867307f4d80f 1108:5ec23b9bf360
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] = obj.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, 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);
90 s = obj.getBoundarySign(boundary); 93 s = obj.getBoundarySign(boundary);
91 gamm = obj.gamm; 94 gamm = obj.gamm;
92 delt = obj.delt; 95 delt = obj.delt;
93 96
94 97
123 tau = s*a*d1; 126 tau = s*a*d1;
124 sig = -s*a*e; 127 sig = -s*a*e;
125 128
126 closure = obj.Hi*(tau*d2' + sig*d3'); 129 closure = obj.Hi*(tau*d2' + sig*d3');
127 penalty{1} = -obj.Hi*tau; 130 penalty{1} = -obj.Hi*tau;
128 penalty{1} = -obj.Hi*sig; 131 penalty{2} = -obj.Hi*sig;
129 132
130 case 'e' 133 case 'e'
131 alpha = obj.alpha; 134 alpha = obj.alpha;
132 tuning = 1.1; 135 tuning = 1.1;
133 136
172 end 175 end
173 176
174 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary, type) 177 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary, type)
175 % u denotes the solution in the own domain 178 % u denotes the solution in the own domain
176 % v denotes the solution in the neighbour domain 179 % v denotes the solution in the neighbour domain
177 [e_u, d1_u, d2_u, d3_u] = obj.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, boundary); 180 e_u = obj.getBoundaryOperator('e', boundary);
181 d1_u = obj.getBoundaryOperator('d1', boundary);
182 d2_u = obj.getBoundaryOperator('d2', boundary);
183 d3_u = obj.getBoundaryOperator('d3', boundary);
178 s_u = obj.getBoundarySign(boundary); 184 s_u = obj.getBoundarySign(boundary);
179 185
180 [e_v, d1_v, d2_v, d3_v] = neighbour_scheme.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, neighbour_boundary); 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);
181 s_v = neighbour_scheme.getBoundarySign(neighbour_boundary); 190 s_v = neighbour_scheme.getBoundarySign(neighbour_boundary);
182 191
183 alpha_u = obj.alpha; 192 alpha_u = obj.alpha;
184 alpha_v = neighbour_scheme.alpha; 193 alpha_v = neighbour_scheme.alpha;
185 194
235 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');
236 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');
237 end 246 end
238 247
239 % Returns the boundary operator op for the boundary specified by the string boundary. 248 % Returns the boundary operator op for the boundary specified by the string boundary.
240 % op -- string or a cell array of strings 249 % op -- string
241 % boundary -- string 250 % boundary -- string
242 function varargout = getBoundaryOperator(obj, op, boundary) 251 function o = getBoundaryOperator(obj, op, boundary)
243 252 assertIsMember(op, {'e', 'd1', 'd2', 'd3'})
244 if ~ismember(boundary, {'l', 'r'}) 253 assertIsMember(boundary, {'l', 'r'})
245 error('No such boundary: boundary = %s',boundary); 254
246 end 255 o = obj.([op, '_', boundary]);
247 256 end
248 if ~iscell(op) 257
249 op = {op}; 258 % Returns square boundary quadrature matrix, of dimension
250 end 259 % corresponding to the number of boundary points
251 260 %
252 for i = 1:numel(op) 261 % boundary -- string
253 switch op{i} 262 % Note: for 1d diffOps, the boundary quadrature is the scalar 1.
254 case 'e' 263 function H_b = getBoundaryQuadrature(obj, boundary)
255 switch boundary 264 assertIsMember(boundary, {'l', 'r'})
256 case 'l' 265
257 e = obj.e_l; 266 H_b = 1;
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 267 end
294 268
295 % Returns the boundary sign. The right boundary is considered the positive boundary 269 % Returns the boundary sign. The right boundary is considered the positive boundary
296 % boundary -- string 270 % boundary -- string
297 function s = getBoundarySign(obj, boundary) 271 function s = getBoundarySign(obj, boundary)
272 assertIsMember(boundary, {'l', 'r'})
273
298 switch boundary 274 switch boundary
299 case {'r'} 275 case {'r'}
300 s = 1; 276 s = 1;
301 case {'l'} 277 case {'l'}
302 s = -1; 278 s = -1;
303 otherwise
304 error('No such boundary: boundary = %s',boundary);
305 end 279 end
306 end 280 end
307 281
308 function N = size(obj) 282 function N = size(obj)
309 N = obj.grid.N; 283 N = obj.grid.N;