Mercurial > repos > public > sbplib

comparison +scheme/Beam2d.m @ 997:78db023a7fe3 feature/getBoundaryOp

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
author Martin Almquist <malmquist@stanford.edu>
date Sat, 12 Jan 2019 11:57:50 -0800
parents a35ed1d124d3
children
comparison
equal deleted inserted replaced
982:a0b3161e44f3 997:78db023a7fe3
119 % neighbour_boundary is a string specifying which boundary to interface to. 119 % neighbour_boundary is a string specifying which boundary to interface to.
120 function [closure, penalty_e,penalty_d] = boundary_condition(obj,boundary,type,data) 120 function [closure, penalty_e,penalty_d] = boundary_condition(obj,boundary,type,data)
121 default_arg('type','dn'); 121 default_arg('type','dn');
122 default_arg('data',0); 122 default_arg('data',0);
123 123
124 [e,d1,d2,d3,s,gamm,delt,halfnorm_inv] = obj.get_boundary_ops(boundary); 124 [e, d1, d2, d3] = obj.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, boundary);
125 s = obj.getBoundarySign(boundary);
126 [gamm, delt] = obj.getBoundaryBorrowing(boundary);
127 halfnorm_inv = obj.getHalfnormInv(boundary);
125 128
126 switch type 129 switch type
127 % Dirichlet-neumann boundary condition 130 % Dirichlet-neumann boundary condition
128 case {'dn'} 131 case {'dn'}
129 alpha = obj.alpha; 132 alpha = obj.alpha;
162 end 165 end
163 166
164 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary, type) 167 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary, type)
165 % u denotes the solution in the own domain 168 % u denotes the solution in the own domain
166 % v denotes the solution in the neighbour domain 169 % v denotes the solution in the neighbour domain
167 [e_u,d1_u,d2_u,d3_u,s_u,gamm_u,delt_u, halfnorm_inv] = obj.get_boundary_ops(boundary); 170 [e_u, d1_u, d2_u, d3_u] = obj.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, boundary);
168 [e_v,d1_v,d2_v,d3_v,s_v,gamm_v,delt_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 171 s_u = obj.getBoundarySign(boundary);
172 [gamm_u, delt_u] = obj.getBoundaryBorrowing(boundary);
173 halfnorm_inv = obj.getHalfnormInv(boundary);
174
175 [e_v, d1_v, d2_v, d3_v] = neighbour_scheme.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, neighbour_boundary);
176 s_v = neighbour_scheme.getBoundarySign(neighbour_boundary);
177 [gamm_v, delt_v] = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
169 178
170 tuning = 2; 179 tuning = 2;
171 180
172 alpha_u = obj.alpha; 181 alpha_u = obj.alpha;
173 alpha_v = neighbour_scheme.alpha; 182 alpha_v = neighbour_scheme.alpha;
190 199
191 closure = halfnorm_inv*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u'); 200 closure = halfnorm_inv*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u');
192 penalty = -halfnorm_inv*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v'); 201 penalty = -halfnorm_inv*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v');
193 end 202 end
194 203
195 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 204 % Returns the boundary operator op for the boundary specified by the string boundary.
196 % The right boundary is considered the positive boundary 205 % op -- string or a cell array of strings
197 function [e,d1,d2,d3,s,gamm, delt, halfnorm_inv] = get_boundary_ops(obj,boundary) 206 % boundary -- string
207 function varargout = getBoundaryOperator(obj, op, boundary)
208
209 if ~iscell(op)
210 op = {op};
211 end
212
213 for i = 1:numel(op)
214 switch op{i}
215 case 'e'
216 switch boundary
217 case 'w'
218 e = obj.e_w;
219 case 'e'
220 e = obj.e_e;
221 case 's'
222 e = obj.e_s;
223 case 'n'
224 e = obj.e_n;
225 otherwise
226 error('No such boundary: boundary = %s',boundary);
227 end
228 varargout{i} = e;
229
230 case 'd1'
231 switch boundary
232 case 'w'
233 d1 = obj.d1_w;
234 case 'e'
235 d1 = obj.d1_e;
236 case 's'
237 d1 = obj.d1_s;
238 case 'n'
239 d1 = obj.d1_n;
240 otherwise
241 error('No such boundary: boundary = %s',boundary);
242 end
243 varargout{i} = d1;
244 end
245
246 case 'd2'
247 switch boundary
248 case 'w'
249 d2 = obj.d2_w;
250 case 'e'
251 d2 = obj.d2_e;
252 case 's'
253 d2 = obj.d2_s;
254 case 'n'
255 d2 = obj.d2_n;
256 otherwise
257 error('No such boundary: boundary = %s',boundary);
258 end
259 varargout{i} = d2;
260 end
261
262 case 'd3'
263 switch boundary
264 case 'w'
265 d3 = obj.d3_w;
266 case 'e'
267 d3 = obj.d3_e;
268 case 's'
269 d3 = obj.d3_s;
270 case 'n'
271 d3 = obj.d3_n;
272 otherwise
273 error('No such boundary: boundary = %s',boundary);
274 end
275 varargout{i} = d3;
276 end
277 end
278 end
279
280 % Returns square boundary quadrature matrix, of dimension
281 % corresponding to the number of boundary points
282 %
283 % boundary -- string
284 function H_b = getBoundaryQuadrature(obj, boundary)
285
198 switch boundary 286 switch boundary
199 case 'w' 287 case 'w'
200 e = obj.e_w; 288 H_b = obj.H_y;
201 d1 = obj.d1_w; 289 case 'e'
202 d2 = obj.d2_w; 290 H_b = obj.H_y;
203 d3 = obj.d3_w; 291 case 's'
292 H_b = obj.H_x;
293 case 'n'
294 H_b = obj.H_x;
295 otherwise
296 error('No such boundary: boundary = %s',boundary);
297 end
298 end
299
300 % Returns the boundary sign. The right boundary is considered the positive boundary
301 % boundary -- string
302 function s = getBoundarySign(obj, boundary)
303 switch boundary
304 case {'e','n'}
305 s = 1;
306 case {'w','s'}
204 s = -1; 307 s = -1;
308 otherwise
309 error('No such boundary: boundary = %s',boundary);
310 end
311 end
312
313 % Returns the halfnorm_inv used in SATs. TODO: better notation
314 function Hinv = getHalfnormInv(obj, boundary)
315 switch boundary
316 case 'w'
317 Hinv = obj.Hix;
318 case 'e'
319 Hinv = obj.Hix;
320 case 's'
321 Hinv = obj.Hiy;
322 case 'n'
323 Hinv = obj.Hiy;
324 otherwise
325 error('No such boundary: boundary = %s',boundary);
326 end
327 end
328
329 % Returns borrowing constant gamma
330 % boundary -- string
331 function [gamm, delt] = getBoundaryBorrowing(obj, boundary)
332 switch boundary
333 case {'w','e'}
205 gamm = obj.gamm_x; 334 gamm = obj.gamm_x;
206 delt = obj.delt_x; 335 delt = obj.delt_x;
207 halfnorm_inv = obj.Hix; 336 case {'s','n'}
208 case 'e'
209 e = obj.e_e;
210 d1 = obj.d1_e;
211 d2 = obj.d2_e;
212 d3 = obj.d3_e;
213 s = 1;
214 gamm = obj.gamm_x;
215 delt = obj.delt_x;
216 halfnorm_inv = obj.Hix;
217 case 's'
218 e = obj.e_s;
219 d1 = obj.d1_s;
220 d2 = obj.d2_s;
221 d3 = obj.d3_s;
222 s = -1;
223 gamm = obj.gamm_y; 337 gamm = obj.gamm_y;
224 delt = obj.delt_y; 338 delt = obj.delt_y;
225 halfnorm_inv = obj.Hiy;
226 case 'n'
227 e = obj.e_n;
228 d1 = obj.d1_n;
229 d2 = obj.d2_n;
230 d3 = obj.d3_n;
231 s = 1;
232 gamm = obj.gamm_y;
233 delt = obj.delt_y;
234 halfnorm_inv = obj.Hiy;
235 otherwise 339 otherwise
236 error('No such boundary: boundary = %s',boundary); 340 error('No such boundary: boundary = %s',boundary);
237 end 341 end
238 end 342 end
239 343