Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger.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 0c504a21432d
children
comparison
equal deleted inserted replaced
1196:f6c571d8f22f 1197:433c89bf19e0
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 assertIsMember(boundary, {'l', 'r'})
118
119 if ~iscell(op)
120 op = {op};
121 end
122
123 for i = 1:numel(op)
124 switch op{i}
125 case 'e'
126 switch boundary
127 case 'l'
128 e = obj.e_l;
129 case 'r'
130 e = obj.e_r;
131 end
132 varargout{i} = e;
133
134 case 'd'
135 switch boundary
136 case 'l'
137 d = obj.d1_l;
138 case 'r'
139 d = obj.d1_r;
140 end
141 varargout{i} = d;
142 end
143 end
144 end
145
146 % Returns square boundary quadrature matrix, of dimension
147 % corresponding to the number of boundary points
148 %
149 % boundary -- string
150 % Note: for 1d diffOps, the boundary quadrature is the scalar 1.
151 function H_b = getBoundaryQuadrature(obj, boundary)
152 assertIsMember(boundary, {'l', 'r'})
153
154 H_b = 1;
155 end
156
157 % Returns the boundary sign. The right boundary is considered the positive boundary
158 % boundary -- string
159 function s = getBoundarySign(obj, boundary)
160 assertIsMember(boundary, {'l', 'r'})
161
112 switch boundary 162 switch boundary
113 case 'l' 163 case {'r'}
114 e = obj.e_l; 164 s = 1;
115 d = obj.d1_l; 165 case {'l'}
116 s = -1; 166 s = -1;
117 case 'r'
118 e = obj.e_r;
119 d = obj.d1_r;
120 s = 1;
121 otherwise
122 error('No such boundary: boundary = %s',boundary);
123 end 167 end
124 end 168 end
125 169
126 function N = size(obj) 170 function N = size(obj)
127 N = obj.m; 171 N = obj.m;
128 end 172 end
129 173
130 end 174 end
131
132 methods(Static)
133 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
134 % and bound_v of scheme schm_v.
135 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
136 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
137 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
138 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
139 end
140 end
141 end 175 end