Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2d.m @ 939:46f5dc61d90b feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update Schrodinger2d to work with new interface types, similar to LaplCurv.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 04 Dec 2018 14:25:03 -0800
parents 1cdf5ead2a16
children 31186559236d
comparison
equal deleted inserted replaced
938:97291e1bd57c 939:46f5dc61d90b
196 otherwise 196 otherwise
197 error('No such boundary condition: type = %s',type); 197 error('No such boundary condition: type = %s',type);
198 end 198 end
199 end 199 end
200 200
201 % opts Struct that specifies the interface coupling. 201 % type Struct that specifies the interface coupling.
202 % Fields: 202 % Fields:
203 % -- interpolation: struct of interpolation operators (empty for conforming grids) 203 % -- interpolation: type of interpolation, default 'none'
204 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,opts) 204 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
205 if isempty(opts) 205
206 defaultType.interpolation = 'none';
207 default_struct('type', defaultType);
208
209 switch type.interpolation
210 case {'none', ''}
206 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary); 211 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary);
207 else 212 case {'op','OP'}
208 assertType(opts, 'struct'); 213 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
209 if isfield(opts, 'I_local2neighbor') 214 otherwise
210 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts); 215 error('Unknown type of interpolation: %s ', type.interpolation);
211 else 216 end
212 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,opts); 217 end
213 end 218
214 end 219 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
215 end
216
217 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
218 % u denotes the solution in the own domain 220 % u denotes the solution in the own domain
219 % v denotes the solution in the neighbour domain 221 % v denotes the solution in the neighbour domain
220 222
221 % Get boundary operators 223 % Get boundary operators
222 [e_neighbour, d_neighbour] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 224 [e_neighbour, d_neighbour] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
235 penalty = -tau*Hi*d*H_gamma*e_neighbour' ... 237 penalty = -tau*Hi*d*H_gamma*e_neighbour' ...
236 -sigma*Hi*e*H_gamma*d_neighbour'; 238 -sigma*Hi*e*H_gamma*d_neighbour';
237 239
238 end 240 end
239 241
240 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts) 242 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
243
244 % User can request special interpolation operators by specifying type.interpOpSet
245 default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
246 interpOpSet = type.interpOpSet;
247
241 % u denotes the solution in the own domain 248 % u denotes the solution in the own domain
242 % v denotes the solution in the neighbour domain 249 % v denotes the solution in the neighbour domain
243 250 [e_v, d_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
244 % Get boundary operators 251 [e_u, d_u, H_gamma] = obj.get_boundary_ops(boundary);
245 [e_neighbour, d_neighbour] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
246 [e, d, H_gamma] = obj.get_boundary_ops(boundary);
247 Hi = obj.Hi; 252 Hi = obj.Hi;
248 a = obj.a; 253 a = obj.a;
249 254
250 % Get outward unit normal component 255 % Get outward unit normal component
251 [~, n] = obj.get_boundary_number(boundary); 256 [~, n] = obj.get_boundary_number(boundary);
252 257
253 % Get interpolation operators from opts 258
254 I_u2v_good = opts.I_local2neighbor.good; 259 % Find the number of grid points along the interface
255 I_u2v_bad = opts.I_local2neighbor.bad; 260 m_u = size(e_u, 2);
256 I_v2u_good = opts.I_neighbor2local.good; 261 m_v = size(e_v, 2);
257 I_v2u_bad = opts.I_neighbor2local.bad; 262
263 % Build interpolation operators
264 intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order);
265 Iu2v = intOps.Iu2v;
266 Iv2u = intOps.Iv2u;
258 267
259 sigma = -n*1i*a/2; 268 sigma = -n*1i*a/2;
260 tau = -n*(1i*a)'/2; 269 tau = -n*(1i*a)'/2;
261 270
262 closure = tau*Hi*d*H_gamma*e' + sigma*Hi*e*H_gamma*d'; 271 closure = tau*Hi*d_u*H_gamma*e_u' + sigma*Hi*e_u*H_gamma*d_u';
263 penalty = -tau*Hi*d*H_gamma*I_v2u_good*e_neighbour' ... 272 penalty = -tau*Hi*d_u*H_gamma*Iv2u.good*e_v' ...
264 -sigma*Hi*e*H_gamma*I_v2u_bad*d_neighbour'; 273 -sigma*Hi*e_u*H_gamma*Iv2u.bad*d_v';
265 274
266 end 275 end
267 276
268 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 277 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
269 function [j, nj] = get_boundary_number(obj, boundary) 278 function [j, nj] = get_boundary_number(obj, boundary)