Mercurial > repos > public > sbplib

comparison +scheme/Utux2d.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 8a9393084b30 84200bbae101
children
comparison
equal deleted inserted replaced
1196:f6c571d8f22f 1197:433c89bf19e0
9 % Can either be a constant vector or a cell array of function handles. 9 % Can either be a constant vector or a cell array of function handles.
10 10
11 H % Discrete norm 11 H % Discrete norm
12 H_x, H_y % Norms in the x and y directions 12 H_x, H_y % Norms in the x and y directions
13 Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms 13 Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms
14 H_w, H_e, H_s, H_n % Boundary quadratures
14 15
15 % Derivatives 16 % Derivatives
16 Dx, Dy 17 Dx, Dy
17 18
18 % Boundary operators 19 % Boundary operators
57 Hx = ops_x.H; 58 Hx = ops_x.H;
58 Hy = ops_y.H; 59 Hy = ops_y.H;
59 Hxi = ops_x.HI; 60 Hxi = ops_x.HI;
60 Hyi = ops_y.HI; 61 Hyi = ops_y.HI;
61 62
63 obj.H_w = Hy;
64 obj.H_e = Hy;
65 obj.H_s = Hx;
66 obj.H_n = Hx;
62 obj.H_x = Hx; 67 obj.H_x = Hx;
63 obj.H_y = Hy; 68 obj.H_y = Hy;
64 obj.H = kron(Hx,Hy); 69 obj.H = kron(Hx,Hy);
65 obj.Hi = kron(Hxi,Hyi); 70 obj.Hi = kron(Hxi,Hyi);
66 obj.Hx = kron(Hx,Iy); 71 obj.Hx = kron(Hx,Iy);
107 % data is a function returning the data that should be applied at the boundary. 112 % data is a function returning the data that should be applied at the boundary.
108 % neighbour_scheme is an instance of Scheme that should be interfaced to. 113 % neighbour_scheme is an instance of Scheme that should be interfaced to.
109 % neighbour_boundary is a string specifying which boundary to interface to. 114 % neighbour_boundary is a string specifying which boundary to interface to.
110 function [closure, penalty] = boundary_condition(obj,boundary,type) 115 function [closure, penalty] = boundary_condition(obj,boundary,type)
111 default_arg('type','dirichlet'); 116 default_arg('type','dirichlet');
112 sigma_left = -1; % Scalar penalty parameter for left boundaries (West/South) 117 s = obj.getBoundarySign(boundary);
113 sigma_right = 1; % Scalar penalty parameter for right boundaries (East/North) 118 e = obj.getBoundaryOperator('e', boundary);
119 H_1d = obj.getOneDirectionalNorm(boundary);
114 switch boundary 120 switch boundary
115 % Can only specify boundary condition where there is inflow 121 % Can only specify boundary condition where there is inflow
116 % Extract the postivie resp. negative part of a, for the left 122 % Extract the postivie resp. negative part of a, for the left
117 % resp. right boundaries, and set other values of a to zero. 123 % resp. right boundaries, and set other values of a to zero.
118 % Then the closure will effectively only contribute to inflow boundaries 124 % Then the closure will effectively only contribute to inflow boundaries
119 case {'w','W','west','West'} 125 case {'w','W','west','West'}
120 a_inflow = obj.a{1}; 126 a_inflow = obj.a{1};
121 a_inflow(a_inflow < 0) = 0; 127 a_inflow(a_inflow < 0) = 0;
122 tau = sigma_left*a_inflow*obj.e_w*obj.H_y;
123 closure = obj.Hi*tau*obj.e_w';
124 case {'e','E','east','East'} 128 case {'e','E','east','East'}
125 a_inflow = obj.a{1}; 129 a_inflow = obj.a{1};
126 a_inflow(a_inflow > 0) = 0; 130 a_inflow(a_inflow > 0) = 0;
127 tau = sigma_right*a_inflow*obj.e_e*obj.H_y;
128 closure = obj.Hi*tau*obj.e_e';
129 case {'s','S','south','South'} 131 case {'s','S','south','South'}
130 a_inflow = obj.a{2}; 132 a_inflow = obj.a{2};
131 a_inflow(a_inflow < 0) = 0; 133 a_inflow(a_inflow < 0) = 0;
132 tau = sigma_left*a_inflow*obj.e_s*obj.H_x;
133 closure = obj.Hi*tau*obj.e_s';
134 case {'n','N','north','North'} 134 case {'n','N','north','North'}
135 a_inflow = obj.a{2}; 135 a_inflow = obj.a{2};
136 a_inflow(a_inflow > 0) = 0; 136 a_inflow(a_inflow > 0) = 0;
137 tau = sigma_right*a_inflow*obj.e_n*obj.H_x; 137 end
138 closure = obj.Hi*tau*obj.e_n'; 138 tau = s*a_inflow*e*H_1d;
139 end 139 closure = obj.Hi*tau*e';
140 penalty = -obj.Hi*tau; 140 penalty = -obj.Hi*tau;
141 end 141 end
142 142
143 % type Struct that specifies the interface coupling. 143 % type Struct that specifies the interface coupling.
144 % Fields: 144 % Fields:
166 166
167 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) 167 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
168 couplingType = type.couplingType; 168 couplingType = type.couplingType;
169 169
170 % Get neighbour boundary operator 170 % Get neighbour boundary operator
171 switch neighbour_boundary 171 e_neighbour = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
172 case {'e','E','east','East'}
173 e_neighbour = neighbour_scheme.e_e;
174 case {'w','W','west','West'}
175 e_neighbour = neighbour_scheme.e_w;
176 case {'n','N','north','North'}
177 e_neighbour = neighbour_scheme.e_n;
178 case {'s','S','south','South'}
179 e_neighbour = neighbour_scheme.e_s;
180 end
181 172
182 switch couplingType 173 switch couplingType
183 174
184 % Upwind coupling (energy dissipation) 175 % Upwind coupling (energy dissipation)
185 case 'upwind' 176 case 'upwind'
224 interpOpSet = type.interpOpSet; 215 interpOpSet = type.interpOpSet;
225 couplingType = type.couplingType; 216 couplingType = type.couplingType;
226 interpolationDamping = type.interpolationDamping; 217 interpolationDamping = type.interpolationDamping;
227 218
228 % Get neighbour boundary operator 219 % Get neighbour boundary operator
229 switch neighbour_boundary 220 e_neighbour = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
230 case {'e','E','east','East'}
231 e_neighbour = neighbour_scheme.e_e;
232 case {'w','W','west','West'}
233 e_neighbour = neighbour_scheme.e_w;
234 case {'n','N','north','North'}
235 e_neighbour = neighbour_scheme.e_n;
236 case {'s','S','south','South'}
237 e_neighbour = neighbour_scheme.e_s;
238 end
239 221
240 switch couplingType 222 switch couplingType
241 223
242 % Upwind coupling (energy dissipation) 224 % Upwind coupling (energy dissipation)
243 case 'upwind' 225 case 'upwind'
317 end 299 end
318 300
319 301
320 end 302 end
321 303
304 % Returns the boundary sign. The right boundary is considered the positive boundary
305 % boundary -- string
306 function s = getBoundarySign(obj, boundary)
307 assertIsMember(boundary, {'w', 'e', 's', 'n'})
308 switch boundary
309 case {'e','n'}
310 s = 1;
311 case {'w','s'}
312 s = -1;
313 end
314 end
315
316 % Returns the boundary operator op for the boundary specified by the string boundary.
317 % op -- string
318 % boundary -- string
319 function o = getBoundaryOperator(obj, op, boundary)
320 assertIsMember(op, {'e'})
321 assertIsMember(boundary, {'w', 'e', 's', 'n'})
322
323 o = obj.([op, '_', boundary]);
324 end
325
326 % Returns square boundary quadrature matrix, of dimension
327 % corresponding to the number of boundary points
328 %
329 % boundary -- string
330 function H_b = getBoundaryQuadrature(obj, boundary)
331 assertIsMember(boundary, {'w', 'e', 's', 'n'})
332
333 H_b = obj.(['H_', boundary]);
334 end
335
336 % Returns square boundary quadrature matrix, of dimension
337 % corresponding to the number of boundary points
338 %
339 % boundary -- string
340 function H_1d = getOneDirectionalNorm(obj, boundary)
341 assertIsMember(boundary, {'w', 'e', 's', 'n'})
342 switch boundary
343 case {'w','e'}
344 H_1d = obj.H_y;
345 case {'s','n'}
346 H_1d = obj.H_x;
347 end
348 end
349
322 function N = size(obj) 350 function N = size(obj)
323 N = obj.m; 351 N = obj.m;
324 end 352 end
325 353
326 end 354 end
327
328 methods(Static)
329 % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u
330 % and bound_v of scheme schm_v.
331 % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l')
332 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
333 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
334 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
335 end
336 end
337 end 355 end