Mercurial > repos > public > sbplib

comparison +scheme/Wave2d.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 5afc774fb7c4
children
comparison
equal deleted inserted replaced
1196:f6c571d8f22f 1197:433c89bf19e0
104 % neighbour_boundary is a string specifying which boundary to interface to. 104 % neighbour_boundary is a string specifying which boundary to interface to.
105 function [closure, penalty] = boundary_condition(obj,boundary,type,data) 105 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
106 default_arg('type','neumann'); 106 default_arg('type','neumann');
107 default_arg('data',0); 107 default_arg('data',0);
108 108
109 [e,d,s,gamm,halfnorm_inv] = obj.get_boundary_ops(boundary); 109 [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary);
110 gamm = obj.getBoundaryBorrowing(boundary);
111 s = obj.getBoundarySign(boundary);
112 halfnorm_inv = obj.getHalfnormInv(boundary);
110 113
111 switch type 114 switch type
112 % Dirichlet boundary condition 115 % Dirichlet boundary condition
113 case {'D','d','dirichlet'} 116 case {'D','d','dirichlet'}
114 alpha = obj.alpha; 117 alpha = obj.alpha;
162 % u denotes the solution in the own domain 165 % u denotes the solution in the own domain
163 % v denotes the solution in the neighbour domain 166 % v denotes the solution in the neighbour domain
164 [e_u,d_u,s_u,gamm_u, halfnorm_inv] = obj.get_boundary_ops(boundary); 167 [e_u,d_u,s_u,gamm_u, halfnorm_inv] = obj.get_boundary_ops(boundary);
165 [e_v,d_v,s_v,gamm_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 168 [e_v,d_v,s_v,gamm_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
166 169
170 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary);
171 gamm_u = obj.getBoundaryBorrowing(boundary);
172 s_u = obj.getBoundarySign(boundary);
173 halfnorm_inv = obj.getHalfnormInv(boundary);
174
175 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary);
176 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
177 s_v = neighbour_scheme.getBoundarySign(neighbour_boundary);
178
167 tuning = 1.1; 179 tuning = 1.1;
168 180
169 alpha_u = obj.alpha; 181 alpha_u = obj.alpha;
170 alpha_v = neighbour_scheme.alpha; 182 alpha_v = neighbour_scheme.alpha;
171 183
181 193
182 closure = halfnorm_inv*( tau*e_u' + sig*alpha_u*d_u'); 194 closure = halfnorm_inv*( tau*e_u' + sig*alpha_u*d_u');
183 penalty = halfnorm_inv*(-tau*e_v' - sig*alpha_v*d_v'); 195 penalty = halfnorm_inv*(-tau*e_v' - sig*alpha_v*d_v');
184 end 196 end
185 197
186 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 198
187 % The right boundary is considered the positive boundary 199 % Returns the boundary operator op for the boundary specified by the string boundary.
188 function [e,d,s,gamm, halfnorm_inv] = get_boundary_ops(obj,boundary) 200 % op -- string or a cell array of strings
201 % boundary -- string
202 function varargout = getBoundaryOperator(obj, op, boundary)
203 assertIsMember(boundary, {'w', 'e', 's', 'n'})
204
205 if ~iscell(op)
206 op = {op};
207 end
208
209 for i = 1:numel(op)
210 switch op{i}
211 case 'e'
212 switch boundary
213 case 'w'
214 e = obj.e_w;
215 case 'e'
216 e = obj.e_e;
217 case 's'
218 e = obj.e_s;
219 case 'n'
220 e = obj.e_n;
221 end
222 varargout{i} = e;
223
224 case 'd'
225 switch boundary
226 case 'w'
227 d = obj.d1_w;
228 case 'e'
229 d = obj.d1_e;
230 case 's'
231 d = obj.d1_s;
232 case 'n'
233 d = obj.d1_n;
234 end
235 varargout{i} = d;
236 end
237 end
238
239 end
240
241 % Returns square boundary quadrature matrix, of dimension
242 % corresponding to the number of boundary points
243 %
244 % boundary -- string
245 function H_b = getBoundaryQuadrature(obj, boundary)
246 assertIsMember(boundary, {'w', 'e', 's', 'n'})
247
189 switch boundary 248 switch boundary
190 case 'w' 249 case 'w'
191 e = obj.e_w; 250 H_b = obj.H_y;
192 d = obj.d1_w; 251 case 'e'
252 H_b = obj.H_y;
253 case 's'
254 H_b = obj.H_x;
255 case 'n'
256 H_b = obj.H_x;
257 end
258 end
259
260 % Returns borrowing constant gamma
261 % boundary -- string
262 function gamm = getBoundaryBorrowing(obj, boundary)
263 assertIsMember(boundary, {'w', 'e', 's', 'n'})
264
265 switch boundary
266 case {'w','e'}
267 gamm = obj.gamm_x;
268 case {'s','n'}
269 gamm = obj.gamm_y;
270 end
271 end
272
273 % Returns the boundary sign. The right boundary is considered the positive boundary
274 % boundary -- string
275 function s = getBoundarySign(obj, boundary)
276 assertIsMember(boundary, {'w', 'e', 's', 'n'})
277
278 switch boundary
279 case {'e','n'}
280 s = 1;
281 case {'w','s'}
193 s = -1; 282 s = -1;
194 gamm = obj.gamm_x; 283 end
195 halfnorm_inv = obj.Hix; 284 end
285
286 % Returns the halfnorm_inv used in SATs. TODO: better notation
287 function Hinv = getHalfnormInv(obj, boundary)
288 assertIsMember(boundary, {'w', 'e', 's', 'n'})
289
290 switch boundary
291 case 'w'
292 Hinv = obj.Hix;
196 case 'e' 293 case 'e'
197 e = obj.e_e; 294 Hinv = obj.Hix;
198 d = obj.d1_e;
199 s = 1;
200 gamm = obj.gamm_x;
201 halfnorm_inv = obj.Hix;
202 case 's' 295 case 's'
203 e = obj.e_s; 296 Hinv = obj.Hiy;
204 d = obj.d1_s;
205 s = -1;
206 gamm = obj.gamm_y;
207 halfnorm_inv = obj.Hiy;
208 case 'n' 297 case 'n'
209 e = obj.e_n; 298 Hinv = obj.Hiy;
210 d = obj.d1_n;
211 s = 1;
212 gamm = obj.gamm_y;
213 halfnorm_inv = obj.Hiy;
214 otherwise
215 error('No such boundary: boundary = %s',boundary);
216 end 299 end
217 end 300 end
218 301
219 function N = size(obj) 302 function N = size(obj)
220 N = prod(obj.m); 303 N = prod(obj.m);
221 end 304 end
222 305
223 end 306 end
224
225 methods(Static)
226 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
227 % and bound_v of scheme schm_v.
228 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
229 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
230 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
231 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
232 end
233 end
234 end 307 end