Mercurial > repos > public > sbplib
changeset 1022:234c1c02ea39 feature/advectionRV
Add dirichlet bc for north and south boundary and handle cases where the wave speed changes in sign.
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Mon, 07 Jan 2019 12:06:06 +0100 |
| parents | cc61dde120cd |
| children | defc9d0cc1f2 |
| files | +scheme/Utux2d.m |
| diffstat | 1 files changed, 27 insertions(+), 8 deletions(-) [+] |
line wrap: on
line diff
diff -r cc61dde120cd -r 234c1c02ea39 +scheme/Utux2d.m --- a/+scheme/Utux2d.m Wed Dec 19 20:00:27 2018 +0100 +++ b/+scheme/Utux2d.m Mon Jan 07 12:06:06 2019 +0100 @@ -105,30 +105,49 @@ else obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy); end - end % Closure functions return the opertors applied to the own domain to close the boundary % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. - % type is a string specifying the type of boundary condition if there are several. + % type is a string specifying the type of boundary condition if there are several. %TBD Remove type here? Only dirichlet applicable? % data is a function returning the data that should be applied at the boundary. % neighbour_scheme is an instance of Scheme that should be interfaced to. % neighbour_boundary is a string specifying which boundary to interface to. function [closure, penalty] = boundary_condition(obj,boundary,type) default_arg('type','dirichlet'); - - sigma = -1; % Scalar penalty parameter + sigma_left = -1; % Scalar penalty parameter for left boundaries + sigma_right = 1; % Scalar penalty parameter for right boundaries switch boundary case {'w','W','west','West'} - tau = sigma*obj.a{1}*obj.e_w*obj.H_y; + % Can only specify boundary condition where there is inflow + % Extract positive part of a{1} + a_pos = obj.a{1}; + a_pos(a_pos < 0) = 0; + tau = sigma_left*a_pos*obj.e_w*obj.H_y; closure = obj.Hi*tau*obj.e_w'; - case {'s','S','south','South'} - tau = sigma*obj.a{2}*obj.e_s*obj.H_x; + % Can only specify boundary condition where there is inflow + % Extract positive part of a{2} + a_pos = obj.a{2}; + a_pos(a_pos < 0) = 0; + tau = sigma_left*a_pos*obj.e_s*obj.H_x; closure = obj.Hi*tau*obj.e_s'; + case {'e','E','east','East'} + % Can only specify boundary condition where there is inflow + % Extract negative part of a{1} + a_neg = obj.a{1}; + a_neg(a_neg > 0) = 0; + tau = sigma_right*a_neg*obj.e_e*obj.H_y; + closure = obj.Hi*tau*obj.e_e'; + case {'n','N','north','North'} + % Can only specify boundary condition where there is inflow + % Extract negative part of a{2} + a_neg = obj.a{2}; + a_neg(a_neg > 0) = 0; + tau = sigma_right*a_neg*obj.e_n*obj.H_x; + closure = obj.Hi*tau*obj.e_n'; end penalty = -obj.Hi*tau; - end % type Struct that specifies the interface coupling.