Mercurial > repos > public > sbplib

comparison +scheme/Wave2dCurve.m @ 347:85c2fe06d551 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implemented characteristic form BC in Wave2dCurve.
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 16 Nov 2016 16:16:44 -0800
parents e070ebd94d9d
children 447ceb41fb65
comparison
equal deleted inserted replaced
346:33b2ef863519 347:85c2fe06d551
153 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 153 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
154 % type is a string specifying the type of boundary condition if there are several. 154 % type is a string specifying the type of boundary condition if there are several.
155 % data is a function returning the data that should be applied at the boundary. 155 % data is a function returning the data that should be applied at the boundary.
156 % neighbour_scheme is an instance of Scheme that should be interfaced to. 156 % neighbour_scheme is an instance of Scheme that should be interfaced to.
157 % neighbour_boundary is a string specifying which boundary to interface to. 157 % neighbour_boundary is a string specifying which boundary to interface to.
158 function [closure, penalty] = boundary_condition(obj,boundary,type) 158 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
159 default_arg('type','neumann'); 159 default_arg('type','neumann');
160 default_arg('parameter', []);
160 161
161 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv] = obj.get_boundary_ops(boundary); 162 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv] = obj.get_boundary_ops(boundary);
162 163
163 switch type 164 switch type
164 % Dirichlet boundary condition 165 % Dirichlet boundary condition
199 tau1 = -s; 200 tau1 = -s;
200 tau2 = 0; 201 tau2 = 0;
201 tau = c.^2 * obj.Ji*(tau1*e + tau2*d); 202 tau = c.^2 * obj.Ji*(tau1*e + tau2*d);
202 203
203 closure = halfnorm_inv*tau*d'; 204 closure = halfnorm_inv*tau*d';
204 penalty = halfnorm_inv*tau; 205 penalty = -halfnorm_inv*tau;
205 206
207 % Characteristic boundary condition
208 case {'characteristic', 'char', 'c'}
209 default_arg('parameter', 1);
210 beta = parameter;
211 c = obj.c;
212
213 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));
215 d = s*(a_n * d_n' + a_t*d_t')'; % outward facing normal derivative
216
217 tau = -c.^2 * 1/beta * obj.Ji*e;
218
219 closure{1} = halfnorm_inv*tau*e';
220 closure{2} = halfnorm_inv*tau*beta*d';
221 penalty = -halfnorm_inv*tau;
206 222
207 % Unknown, boundary condition 223 % Unknown, boundary condition
208 otherwise 224 otherwise
209 error('No such boundary condition: type = %s',type); 225 error('No such boundary condition: type = %s',type);
210 end 226 end