Mercurial > repos > public > sbplib

comparison +scheme/Wave2d.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 706d1c2b4199
children 8d73fcdb07a5
comparison
equal deleted inserted replaced
982:a0b3161e44f3 997:78db023a7fe3
104 % neighbour_boundary is a string specifying which boundary to interface to. 104 % neighbour_boundary is a string specifying which boundary to interface to.
105 function [closure, penalty] = boundary_condition(obj,boundary,type,data) 105 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
106 default_arg('type','neumann'); 106 default_arg('type','neumann');
107 default_arg('data',0); 107 default_arg('data',0);
108 108
109 [e,d,s,gamm,halfnorm_inv] = obj.get_boundary_ops(boundary); 109 [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary);
110 gamm = obj.getBoundaryBorrowing(boundary);
111 s = obj.getBoundarySign(boundary);
112 halfnorm_inv = obj.getHalfnormInv(boundary);
110 113
111 switch type 114 switch type
112 % Dirichlet boundary condition 115 % Dirichlet boundary condition
113 case {'D','d','dirichlet'} 116 case {'D','d','dirichlet'}
114 alpha = obj.alpha; 117 alpha = obj.alpha;
162 % u denotes the solution in the own domain 165 % u denotes the solution in the own domain
163 % v denotes the solution in the neighbour domain 166 % v denotes the solution in the neighbour domain
164 [e_u,d_u,s_u,gamm_u, halfnorm_inv] = obj.get_boundary_ops(boundary); 167 [e_u,d_u,s_u,gamm_u, halfnorm_inv] = obj.get_boundary_ops(boundary);
165 [e_v,d_v,s_v,gamm_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 168 [e_v,d_v,s_v,gamm_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
166 169
170 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary);
171 gamm_u = obj.getBoundaryBorrowing(boundary);
172 s_u = obj.getBoundarySign(boundary);
173 halfnorm_inv = obj.getHalfnormInv(boundary);
174
175 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary);
176 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
177 s_v = neighbour_scheme.getBoundarySign(neighbour_boundary);
178
167 tuning = 1.1; 179 tuning = 1.1;
168 180
169 alpha_u = obj.alpha; 181 alpha_u = obj.alpha;
170 alpha_v = neighbour_scheme.alpha; 182 alpha_v = neighbour_scheme.alpha;
171 183
181 193
182 closure = halfnorm_inv*( tau*e_u' + sig*alpha_u*d_u'); 194 closure = halfnorm_inv*( tau*e_u' + sig*alpha_u*d_u');
183 penalty = halfnorm_inv*(-tau*e_v' - sig*alpha_v*d_v'); 195 penalty = halfnorm_inv*(-tau*e_v' - sig*alpha_v*d_v');
184 end 196 end
185 197
186 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 198
187 % The right boundary is considered the positive boundary 199 % Returns the boundary operator op for the boundary specified by the string boundary.
188 function [e,d,s,gamm, halfnorm_inv] = get_boundary_ops(obj,boundary) 200 % op -- string or a cell array of strings
201 % boundary -- string
202 function varargout = getBoundaryOperator(obj, op, boundary)
203
204 if ~iscell(op)
205 op = {op};
206 end
207
208 for i = 1:numel(op)
209 switch op{i}
210 case 'e'
211 switch boundary
212 case 'w'
213 e = obj.e_w;
214 case 'e'
215 e = obj.e_e;
216 case 's'
217 e = obj.e_s;
218 case 'n'
219 e = obj.e_n;
220 otherwise
221 error('No such boundary: boundary = %s',boundary);
222 end
223 varargout{i} = e;
224
225 case 'd'
226 switch boundary
227 case 'w'
228 d = obj.d1_w;
229 case 'e'
230 d = obj.d1_e;
231 case 's'
232 d = obj.d1_s;
233 case 'n'
234 d = obj.d1_n;
235 otherwise
236 error('No such boundary: boundary = %s',boundary);
237 end
238 varargout{i} = d;
239 end
240 end
241
242 end
243
244 % Returns square boundary quadrature matrix, of dimension
245 % corresponding to the number of boundary points
246 %
247 % boundary -- string
248 function H_b = getBoundaryQuadrature(obj, boundary)
249
189 switch boundary 250 switch boundary
190 case 'w' 251 case 'w'
191 e = obj.e_w; 252 H_b = obj.H_y;
192 d = obj.d1_w; 253 case 'e'
254 H_b = obj.H_y;
255 case 's'
256 H_b = obj.H_x;
257 case 'n'
258 H_b = obj.H_x;
259 otherwise
260 error('No such boundary: boundary = %s',boundary);
261 end
262 end
263
264 % Returns borrowing constant gamma
265 % boundary -- string
266 function gamm = getBoundaryBorrowing(obj, boundary)
267 switch boundary
268 case {'w','e'}
269 gamm = obj.gamm_x;
270 case {'s','n'}
271 gamm = obj.gamm_y;
272 otherwise
273 error('No such boundary: boundary = %s',boundary);
274 end
275 end
276
277 % Returns the boundary sign. The right boundary is considered the positive boundary
278 % boundary -- string
279 function s = getBoundarySign(obj, boundary)
280 switch boundary
281 case {'e','n'}
282 s = 1;
283 case {'w','s'}
193 s = -1; 284 s = -1;
194 gamm = obj.gamm_x; 285 otherwise
195 halfnorm_inv = obj.Hix; 286 error('No such boundary: boundary = %s',boundary);
287 end
288 end
289
290 % Returns the halfnorm_inv used in SATs. TODO: better notation
291 function Hinv = getHalfnormInv(obj, boundary)
292 switch boundary
293 case 'w'
294 Hinv = obj.Hix;
196 case 'e' 295 case 'e'
197 e = obj.e_e; 296 Hinv = obj.Hix;
198 d = obj.d1_e;
199 s = 1;
200 gamm = obj.gamm_x;
201 halfnorm_inv = obj.Hix;
202 case 's' 297 case 's'
203 e = obj.e_s; 298 Hinv = obj.Hiy;
204 d = obj.d1_s;
205 s = -1;
206 gamm = obj.gamm_y;
207 halfnorm_inv = obj.Hiy;
208 case 'n' 299 case 'n'
209 e = obj.e_n; 300 Hinv = obj.Hiy;
210 d = obj.d1_n;
211 s = 1;
212 gamm = obj.gamm_y;
213 halfnorm_inv = obj.Hiy;
214 otherwise 301 otherwise
215 error('No such boundary: boundary = %s',boundary); 302 error('No such boundary: boundary = %s',boundary);
216 end 303 end
217 end 304 end
218 305