Mercurial > repos > public > sbplib

comparison +scheme/Utux.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 0c504a21432d
children
comparison
equal deleted inserted replaced
1196:f6c571d8f22f 1197:433c89bf19e0
66 % data is a function returning the data that should be applied at the boundary. 66 % data is a function returning the data that should be applied at the boundary.
67 % neighbour_scheme is an instance of Scheme that should be interfaced to. 67 % neighbour_scheme is an instance of Scheme that should be interfaced to.
68 % neighbour_boundary is a string specifying which boundary to interface to. 68 % neighbour_boundary is a string specifying which boundary to interface to.
69 function [closure, penalty] = boundary_condition(obj,boundary,type) 69 function [closure, penalty] = boundary_condition(obj,boundary,type)
70 default_arg('type','dirichlet'); 70 default_arg('type','dirichlet');
71 sigma_left = -1; % Scalar penalty parameter for left boundary 71 s = obj.getBoundarySign(boundary);
72 sigma_right = 1; % Scalar penalty parameter for right boundary 72 e = obj.getBoundaryOperator('e', boundary);
73 switch boundary 73 switch boundary
74 % Can only specify boundary condition where there is inflow 74 % Can only specify boundary condition where there is inflow
75 % Extract the postivie resp. negative part of a, for the left 75 % Extract the postivie resp. negative part of a, for the left
76 % resp. right boundary, and set other values of a to zero. 76 % resp. right boundary, and set other values of a to zero.
77 % Then the closure will effectively only contribute to inflow boundaries 77 % Then the closure will effectively only contribute to inflow boundaries
78 case {'l','L','left','Left'} 78 case {'l'}
79 a_inflow = obj.a; 79 a_inflow = obj.a;
80 a_inflow(a_inflow < 0) = 0; 80 a_inflow(a_inflow < 0) = 0;
81 tau = sigma_left*a_inflow*obj.e_l; 81 case {'r'}
82 closure = obj.Hi*tau*obj.e_l';
83 case {'r','R','right','Right'}
84 a_inflow = obj.a; 82 a_inflow = obj.a;
85 a_inflow(a_inflow > 0) = 0; 83 a_inflow(a_inflow > 0) = 0;
86 tau = sigma_right*a_inflow*obj.e_r;
87 closure = obj.Hi*tau*obj.e_r';
88 end 84 end
85 tau = s*a_inflow*e;
86 closure = obj.Hi*tau*e';
89 penalty = -obj.Hi*tau; 87 penalty = -obj.Hi*tau;
90
91 end 88 end
92 89
93 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) 90 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
94 switch boundary 91 switch boundary
95 % Upwind coupling 92 % Upwind coupling
103 penalty = -obj.Hi*tau*neighbour_scheme.e_l'; 100 penalty = -obj.Hi*tau*neighbour_scheme.e_l';
104 end 101 end
105 102
106 end 103 end
107 104
105 % Returns the boundary sign. The right boundary is considered the positive boundary
106 % boundary -- string
107 function s = getBoundarySign(obj, boundary)
108 assertIsMember(boundary, {'l', 'r'})
109
110 switch boundary
111 case {'r'}
112 s = 1;
113 case {'l'}
114 s = -1;
115 end
116 end
117
118 % Returns the boundary operator op for the boundary specified by the string boundary.
119 % op -- string
120 % boundary -- string
121 function o = getBoundaryOperator(obj, op, boundary)
122 assertIsMember(op, {'e'})
123 assertIsMember(boundary, {'l', 'r'})
124
125 o = obj.([op, '_', boundary]);
126 end
127
128 % Returns square boundary quadrature matrix, of dimension
129 % corresponding to the number of boundary points
130 %
131 % boundary -- string
132 % Note: for 1d diffOps, the boundary quadrature is the scalar 1.
133 function H_b = getBoundaryQuadrature(obj, boundary)
134 assertIsMember(boundary, {'l', 'r'})
135
136 H_b = 1;
137 end
138
108 function N = size(obj) 139 function N = size(obj)
109 N = obj.m; 140 N = obj.m;
110 end 141 end
111 142
112 end 143 end
113
114 methods(Static)
115 % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u
116 % and bound_v of scheme schm_v.
117 % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l')
118 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
119 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
120 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
121 end
122 end
123 end 144 end