Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2d.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 3dd7f87c9a1b
children 8d73fcdb07a5
comparison
equal deleted inserted replaced
982:a0b3161e44f3 997:78db023a7fe3
160 % neighbour_boundary is a string specifying which boundary to interface to. 160 % neighbour_boundary is a string specifying which boundary to interface to.
161 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 161 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
162 default_arg('type','Neumann'); 162 default_arg('type','Neumann');
163 default_arg('parameter', []); 163 default_arg('parameter', []);
164 164
165 % j is the coordinate direction of the boundary
166 % nj: outward unit normal component. 165 % nj: outward unit normal component.
167 % nj = -1 for west, south, bottom boundaries 166 % nj = -1 for west, south, bottom boundaries
168 % nj = 1 for east, north, top boundaries 167 % nj = 1 for east, north, top boundaries
169 [j, nj] = obj.get_boundary_number(boundary); 168 nj = obj.getBoundarySign(boundary);
170 switch nj 169 [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary);
171 case 1 170 H_gamma = obj.getBoundaryQuadrature(boundary);
172 e = obj.e_r;
173 d = obj.d1_r;
174 case -1
175 e = obj.e_l;
176 d = obj.d1_l;
177 end
178
179 Hi = obj.Hi; 171 Hi = obj.Hi;
180 H_gamma = obj.H_boundary{j}; 172 a = e'*obj.a*e;
181 a = e{j}'*obj.a*e{j};
182 173
183 switch type 174 switch type
184 175
185 % Dirichlet boundary condition 176 % Dirichlet boundary condition
186 case {'D','d','dirichlet','Dirichlet'} 177 case {'D','d','dirichlet','Dirichlet'}
187 closure = nj*Hi*d{j}*a*1i*H_gamma*(e{j}' ); 178 closure = nj*Hi*d*a*1i*H_gamma*(e' );
188 penalty = -nj*Hi*d{j}*a*1i*H_gamma; 179 penalty = -nj*Hi*d*a*1i*H_gamma;
189 180
190 % Free boundary condition 181 % Free boundary condition
191 case {'N','n','neumann','Neumann'} 182 case {'N','n','neumann','Neumann'}
192 closure = -nj*Hi*e{j}*a*1i*H_gamma*(d{j}' ); 183 closure = -nj*Hi*e*a*1i*H_gamma*(d' );
193 penalty = nj*Hi*e{j}*a*1i*H_gamma; 184 penalty = nj*Hi*e*a*1i*H_gamma;
194 185
195 % Unknown boundary condition 186 % Unknown boundary condition
196 otherwise 187 otherwise
197 error('No such boundary condition: type = %s',type); 188 error('No such boundary condition: type = %s',type);
198 end 189 end
219 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) 210 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
220 % u denotes the solution in the own domain 211 % u denotes the solution in the own domain
221 % v denotes the solution in the neighbour domain 212 % v denotes the solution in the neighbour domain
222 213
223 % Get boundary operators 214 % Get boundary operators
224 [e_neighbour, d_neighbour] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 215 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary);
225 [e, d, H_gamma] = obj.get_boundary_ops(boundary); 216 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary);
217 H_gamma = obj.getBoundaryQuadrature(boundary);
226 Hi = obj.Hi; 218 Hi = obj.Hi;
227 a = obj.a; 219 a = obj.a;
228 220
229 % Get outward unit normal component 221 % Get outward unit normal component
230 [~, n] = obj.get_boundary_number(boundary); 222 n = obj.getBoundarySign(boundary);
231 223
232 Hi = obj.Hi; 224 Hi = obj.Hi;
233 sigma = -n*1i*a/2; 225 sigma = -n*1i*a/2;
234 tau = -n*(1i*a)'/2; 226 tau = -n*(1i*a)'/2;
235 227
245 default_field(type, 'interpOpSet', @sbp.InterpOpsOP); 237 default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
246 interpOpSet = type.interpOpSet; 238 interpOpSet = type.interpOpSet;
247 239
248 % u denotes the solution in the own domain 240 % u denotes the solution in the own domain
249 % v denotes the solution in the neighbour domain 241 % v denotes the solution in the neighbour domain
250 [e_v, d_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 242 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary);
251 [e_u, d_u, H_gamma] = obj.get_boundary_ops(boundary); 243 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary);
244 H_gamma = obj.getBoundaryQuadrature(boundary);
252 Hi = obj.Hi; 245 Hi = obj.Hi;
253 a = obj.a; 246 a = obj.a;
254 247
255 % Get outward unit normal component 248 % Get outward unit normal component
256 [~, n] = obj.get_boundary_number(boundary); 249 n = obj.getBoundarySign(boundary);
257 250
258 % Find the number of grid points along the interface 251 % Find the number of grid points along the interface
259 m_u = size(e_u, 2); 252 m_u = size(e_u, 2);
260 m_v = size(e_v, 2); 253 m_v = size(e_v, 2);
261 254
291 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 284 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
292 nj = 1; 285 nj = 1;
293 end 286 end
294 end 287 end
295 288
296 % Returns the boundary ops and sign for the boundary specified by the string boundary. 289 % Returns the boundary operator op for the boundary specified by the string boundary.
297 % The right boundary is considered the positive boundary 290 % op -- string or a cell array of strings
298 function [e, d, H_b] = get_boundary_ops(obj, boundary) 291 % boundary -- string
292 function varargout = getBoundaryOperator(obj, op, boundary)
293
294 if ~iscell(op)
295 op = {op};
296 end
297
298 for i = 1:numel(op)
299 switch op{i}
300 case 'e'
301 switch boundary
302 case 'w'
303 e = obj.e_w;
304 case 'e'
305 e = obj.e_e;
306 case 's'
307 e = obj.e_s;
308 case 'n'
309 e = obj.e_n;
310 otherwise
311 error('No such boundary: boundary = %s',boundary);
312 end
313 varargout{i} = e;
314
315 case 'd'
316 switch boundary
317 case 'w'
318 d = obj.d_w;
319 case 'e'
320 d = obj.d_e;
321 case 's'
322 d = obj.d_s;
323 case 'n'
324 d = obj.d_n;
325 otherwise
326 error('No such boundary: boundary = %s',boundary);
327 end
328 varargout{i} = d;
329 end
330 end
331 end
332
333 % Returns square boundary quadrature matrix, of dimension
334 % corresponding to the number of boundary points
335 %
336 % boundary -- string
337 function H_b = getBoundaryQuadrature(obj, boundary)
299 338
300 switch boundary 339 switch boundary
301 case 'w' 340 case 'w'
302 e = obj.e_w;
303 d = obj.d_w;
304 H_b = obj.H_boundary{1}; 341 H_b = obj.H_boundary{1};
305 case 'e' 342 case 'e'
306 e = obj.e_e;
307 d = obj.d_e;
308 H_b = obj.H_boundary{1}; 343 H_b = obj.H_boundary{1};
309 case 's' 344 case 's'
310 e = obj.e_s;
311 d = obj.d_s;
312 H_b = obj.H_boundary{2}; 345 H_b = obj.H_boundary{2};
313 case 'n' 346 case 'n'
314 e = obj.e_n;
315 d = obj.d_n;
316 H_b = obj.H_boundary{2}; 347 H_b = obj.H_boundary{2};
317 otherwise 348 otherwise
318 error('No such boundary: boundary = %s',boundary); 349 error('No such boundary: boundary = %s',boundary);
319 end 350 end
320 end 351 end
321 352
353 % Returns the boundary sign. The right boundary is considered the positive boundary
354 % boundary -- string
355 function s = getBoundarySign(obj, boundary)
356 switch boundary
357 case {'e','n'}
358 s = 1;
359 case {'w','s'}
360 s = -1;
361 otherwise
362 error('No such boundary: boundary = %s',boundary);
363 end
364 end
365
322 function N = size(obj) 366 function N = size(obj)
323 N = prod(obj.m); 367 N = prod(obj.m);
324 end 368 end
325 end 369 end
326 end 370 end