Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger.m @ 999:337c4d1dcef5 feature/getBoundaryOp

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix overlooked Schrodinger1d scheme.
author Martin Almquist <malmquist@stanford.edu>
date Sat, 12 Jan 2019 13:44:08 -0800
parents 706d1c2b4199
children 8d73fcdb07a5
comparison
equal deleted inserted replaced
998:2b1b944deae1 999:337c4d1dcef5
65 % neighbour_boundary is a string specifying which boundary to interface to. 65 % neighbour_boundary is a string specifying which boundary to interface to.
66 function [closure, penalty] = boundary_condition(obj,boundary,type,data) 66 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
67 default_arg('type','dirichlet'); 67 default_arg('type','dirichlet');
68 default_arg('data',0); 68 default_arg('data',0);
69 69
70 [e,d,s] = obj.get_boundary_ops(boundary); 70 [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary);
71 s = obj.getBoundarySign(boundary);
71 72
72 switch type 73 switch type
73 % Dirichlet boundary condition 74 % Dirichlet boundary condition
74 case {'D','d','dirichlet'} 75 case {'D','d','dirichlet'}
75 tau = s * 1i*d; 76 tau = s * 1i*d;
91 end 92 end
92 93
93 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) 94 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
94 % u denotes the solution in the own domain 95 % u denotes the solution in the own domain
95 % v denotes the solution in the neighbour domain 96 % v denotes the solution in the neighbour domain
96 [e_u,d_u,s_u] = obj.get_boundary_ops(boundary); 97 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary);
97 [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 98 s_u = obj.getBoundarySign(boundary);
99
100 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary);
101 s_v = neighbour_scheme.getBoundarySign(neighbour_boundary);
98 102
99 a = -s_u* 1/2 * 1i ; 103 a = -s_u* 1/2 * 1i ;
100 b = a'; 104 b = a';
101 105
102 tau = b*d_u; 106 tau = b*d_u;
104 108
105 closure = obj.Hi * (tau*e_u' + sig*d_u'); 109 closure = obj.Hi * (tau*e_u' + sig*d_u');
106 penalty = obj.Hi * (-tau*e_v' - sig*d_v'); 110 penalty = obj.Hi * (-tau*e_v' - sig*d_v');
107 end 111 end
108 112
109 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 113 % Returns the boundary operator op for the boundary specified by the string boundary.
110 % The right boundary is considered the positive boundary 114 % op -- string or a cell array of strings
111 function [e,d,s] = get_boundary_ops(obj,boundary) 115 % boundary -- string
116 function varargout = getBoundaryOperator(obj, op, boundary)
117
118 if ~iscell(op)
119 op = {op};
120 end
121
122 for i = 1:numel(op)
123 switch op{i}
124 case 'e'
125 switch boundary
126 case 'l'
127 e = obj.e_l;
128 case 'r'
129 e = obj.e_r;
130 otherwise
131 error('No such boundary: boundary = %s',boundary);
132 end
133 varargout{i} = e;
134
135 case 'd'
136 switch boundary
137 case 'l'
138 d = obj.d1_l;
139 case 'r'
140 d = obj.d1_r;
141 otherwise
142 error('No such boundary: boundary = %s',boundary);
143 end
144 varargout{i} = d;
145 end
146 end
147 end
148
149 % Returns the boundary sign. The right boundary is considered the positive boundary
150 % boundary -- string
151 function s = getBoundarySign(obj, boundary)
112 switch boundary 152 switch boundary
113 case 'l' 153 case {'r'}
114 e = obj.e_l; 154 s = 1;
115 d = obj.d1_l; 155 case {'l'}
116 s = -1; 156 s = -1;
117 case 'r'
118 e = obj.e_r;
119 d = obj.d1_r;
120 s = 1;
121 otherwise 157 otherwise
122 error('No such boundary: boundary = %s',boundary); 158 error('No such boundary: boundary = %s',boundary);
123 end 159 end
124 end 160 end
125 161