Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2d.m @ 1197:433c89bf19e0 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 07 Aug 2019 15:23:42 +0200
parents 8d73fcdb07a5
children
comparison
equal deleted inserted replaced
1196:f6c571d8f22f 1197:433c89bf19e0
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 assertIsMember(boundary, {'w', 'e', 's', 'n'})
294
295 if ~iscell(op)
296 op = {op};
297 end
298
299 for i = 1:numel(op)
300 switch op{i}
301 case 'e'
302 switch boundary
303 case 'w'
304 e = obj.e_w;
305 case 'e'
306 e = obj.e_e;
307 case 's'
308 e = obj.e_s;
309 case 'n'
310 e = obj.e_n;
311 end
312 varargout{i} = e;
313
314 case 'd'
315 switch boundary
316 case 'w'
317 d = obj.d_w;
318 case 'e'
319 d = obj.d_e;
320 case 's'
321 d = obj.d_s;
322 case 'n'
323 d = obj.d_n;
324 end
325 varargout{i} = d;
326 end
327 end
328 end
329
330 % Returns square boundary quadrature matrix, of dimension
331 % corresponding to the number of boundary points
332 %
333 % boundary -- string
334 function H_b = getBoundaryQuadrature(obj, boundary)
335 assertIsMember(boundary, {'w', 'e', 's', 'n'})
299 336
300 switch boundary 337 switch boundary
301 case 'w' 338 case 'w'
302 e = obj.e_w;
303 d = obj.d_w;
304 H_b = obj.H_boundary{1}; 339 H_b = obj.H_boundary{1};
305 case 'e' 340 case 'e'
306 e = obj.e_e;
307 d = obj.d_e;
308 H_b = obj.H_boundary{1}; 341 H_b = obj.H_boundary{1};
309 case 's' 342 case 's'
310 e = obj.e_s;
311 d = obj.d_s;
312 H_b = obj.H_boundary{2}; 343 H_b = obj.H_boundary{2};
313 case 'n' 344 case 'n'
314 e = obj.e_n;
315 d = obj.d_n;
316 H_b = obj.H_boundary{2}; 345 H_b = obj.H_boundary{2};
317 otherwise 346 end
318 error('No such boundary: boundary = %s',boundary); 347 end
348
349 % Returns the boundary sign. The right boundary is considered the positive boundary
350 % boundary -- string
351 function s = getBoundarySign(obj, boundary)
352 assertIsMember(boundary, {'w', 'e', 's', 'n'})
353
354 switch boundary
355 case {'e','n'}
356 s = 1;
357 case {'w','s'}
358 s = -1;
319 end 359 end
320 end 360 end
321 361
322 function N = size(obj) 362 function N = size(obj)
323 N = prod(obj.m); 363 N = prod(obj.m);