Mercurial > repos > public > sbplib

comparison +scheme/Wave2dCurve.m @ 384:32df00102268 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix bug in Characteristic BC for Wave2dCurve.
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 02 Jan 2017 09:24:48 +0100
parents 447ceb41fb65
children 42c89b5eedc0
comparison
equal deleted inserted replaced
383:151b08366d63 384:32df00102268
24 du_e, dv_e 24 du_e, dv_e
25 du_s, dv_s 25 du_s, dv_s
26 du_n, dv_n 26 du_n, dv_n
27 gamm_u, gamm_v 27 gamm_u, gamm_v
28 lambda 28 lambda
29
30 x_u
31 x_v
32 y_u
33 y_v
29 end 34 end
30 35
31 methods 36 methods
32 function obj = Wave2dCurve(g ,order, c, opSet) 37 function obj = Wave2dCurve(g ,order, c, opSet)
33 default_arg('opSet',@sbp.D2Variable); 38 default_arg('opSet',@sbp.D2Variable);
127 obj.du_s = (obj.e_s'*Du)'; 132 obj.du_s = (obj.e_s'*Du)';
128 obj.dv_s = kr(I_u,d1_l_v); 133 obj.dv_s = kr(I_u,d1_l_v);
129 obj.du_n = (obj.e_n'*Du)'; 134 obj.du_n = (obj.e_n'*Du)';
130 obj.dv_n = kr(I_u,d1_r_v); 135 obj.dv_n = kr(I_u,d1_r_v);
131 136
137 obj.x_u = x_u;
138 obj.x_v = x_v;
139 obj.y_u = y_u;
140 obj.y_v = y_v;
141
132 obj.m = m; 142 obj.m = m;
133 obj.h = [h_u h_v]; 143 obj.h = [h_u h_v];
134 obj.order = order; 144 obj.order = order;
135 obj.grid = g; 145 obj.grid = g;
136 146
157 % neighbour_boundary is a string specifying which boundary to interface to. 167 % neighbour_boundary is a string specifying which boundary to interface to.
158 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 168 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
159 default_arg('type','neumann'); 169 default_arg('type','neumann');
160 default_arg('parameter', []); 170 default_arg('parameter', []);
161 171
162 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv] = obj.get_boundary_ops(boundary); 172 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv , ~, ~, ~, scale_factor] = obj.get_boundary_ops(boundary);
163
164 switch type 173 switch type
165 % Dirichlet boundary condition 174 % Dirichlet boundary condition
166 case {'D','d','dirichlet'} 175 case {'D','d','dirichlet'}
167 % v denotes the solution in the neighbour domain 176 % v denotes the solution in the neighbour domain
168 tuning = 1.2; 177 tuning = 1.2;
212 221
213 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); 222 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n));
214 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); 223 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t));
215 d = s*(a_n * d_n' + a_t*d_t')'; % outward facing normal derivative 224 d = s*(a_n * d_n' + a_t*d_t')'; % outward facing normal derivative
216 225
217 tau = -c.^2 * 1/beta * obj.Ji*e; 226 tau = -c.^2 * 1/beta*obj.Ji*e;
218 227
219 closure{1} = halfnorm_inv*tau*e'; 228 closure{1} = halfnorm_inv*tau/c*spdiag(scale_factor)*e';
220 closure{2} = halfnorm_inv*tau*beta*d'; 229 closure{2} = halfnorm_inv*tau*beta*d';
221 penalty = -halfnorm_inv*tau; 230 penalty = -halfnorm_inv*tau;
222 231
223 % Unknown, boundary condition 232 % Unknown, boundary condition
224 otherwise 233 otherwise
265 274
266 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 275 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
267 % The right boundary is considered the positive boundary 276 % The right boundary is considered the positive boundary
268 % 277 %
269 % I -- the indecies of the boundary points in the grid matrix 278 % I -- the indecies of the boundary points in the grid matrix
270 function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I] = get_boundary_ops(obj, boundary) 279 function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor] = get_boundary_ops(obj, boundary)
271 280
272 % gridMatrix = zeros(obj.m(2),obj.m(1)); 281 % gridMatrix = zeros(obj.m(2),obj.m(1));
273 % gridMatrix(:) = 1:numel(gridMatrix); 282 % gridMatrix(:) = 1:numel(gridMatrix);
274 283
275 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); 284 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
282 s = -1; 291 s = -1;
283 292
284 I = ind(1,:); 293 I = ind(1,:);
285 coeff_n = obj.a11(I); 294 coeff_n = obj.a11(I);
286 coeff_t = obj.a12(I); 295 coeff_t = obj.a12(I);
296 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
287 case 'e' 297 case 'e'
288 e = obj.e_e; 298 e = obj.e_e;
289 d_n = obj.du_e; 299 d_n = obj.du_e;
290 d_t = obj.dv_e; 300 d_t = obj.dv_e;
291 s = 1; 301 s = 1;
292 302
293 I = ind(end,:); 303 I = ind(end,:);
294 coeff_n = obj.a11(I); 304 coeff_n = obj.a11(I);
295 coeff_t = obj.a12(I); 305 coeff_t = obj.a12(I);
306 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
296 case 's' 307 case 's'
297 e = obj.e_s; 308 e = obj.e_s;
298 d_n = obj.dv_s; 309 d_n = obj.dv_s;
299 d_t = obj.du_s; 310 d_t = obj.du_s;
300 s = -1; 311 s = -1;
301 312
302 I = ind(:,1)'; 313 I = ind(:,1)';
303 coeff_n = obj.a22(I); 314 coeff_n = obj.a22(I);
304 coeff_t = obj.a12(I); 315 coeff_t = obj.a12(I);
316 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
305 case 'n' 317 case 'n'
306 e = obj.e_n; 318 e = obj.e_n;
307 d_n = obj.dv_n; 319 d_n = obj.dv_n;
308 d_t = obj.du_n; 320 d_t = obj.du_n;
309 s = 1; 321 s = 1;
310 322
311 I = ind(:,end)'; 323 I = ind(:,end)';
312 coeff_n = obj.a22(I); 324 coeff_n = obj.a22(I);
313 coeff_t = obj.a12(I); 325 coeff_t = obj.a12(I);
326 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
314 otherwise 327 otherwise
315 error('No such boundary: boundary = %s',boundary); 328 error('No such boundary: boundary = %s',boundary);
316 end 329 end
317 330
318 switch boundary 331 switch boundary