Mercurial > repos > public > sbplib

comparison +scheme/Burgers2d.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 6cb03209f0a7
children
comparison
equal deleted inserted replaced
1196:f6c571d8f22f 1197:433c89bf19e0
99 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain. 99 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain.
100 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 100 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
101 % type is a string specifying the type of boundary condition if there are several. 101 % type is a string specifying the type of boundary condition if there are several.
102 function [closure, penalty] = boundary_condition(obj,boundary,type) 102 function [closure, penalty] = boundary_condition(obj,boundary,type)
103 default_arg('type','dirichlet'); 103 default_arg('type','dirichlet');
104 [e, H_b, index, s] = obj.get_boundary_ops(boundary); 104 s = obj.getBoundarySign(boundary);
105 e = obj.getBoundaryOperator('e', boundary);
106 indices = obj.getBoundaryIndices(boundary);
107 H_1d = obj.getOneDirectionalNorm(boundary);
105 switch type 108 switch type
106 % Stable dirchlet-like boundary conditions (u+-abs(u))*u/3
107 % with +- at left/right boundaries in each coordinate direction
108 case {'D', 'd', 'dirichlet', 'Dirichlet'} 109 case {'D', 'd', 'dirichlet', 'Dirichlet'}
109 110 penalty_parameter = 1/3;
110 magnitude = 1/3; 111 Tau = s*penalty_parameter*obj.Hi*e*H_1d/2;
111 Tau = s*magnitude*obj.Hi*e*H_b/2; 112 m = obj.grid.m;
112 m = length(index); 113 tau = @(v) Tau*spdiags((v(indices)-s*abs(v(indices))),0,m(1),m(2));
113 tau = @(v) Tau*spdiags((v(index)-s*abs(v(index))),0,m,m); 114 closure = @(v) Tau*((v(indices)-s*abs(v(indices))).*v(indices));
114 closure = @(v) Tau*((v(index)-s*abs(v(index))).*v(index));
115 penalty = @(v) -tau(v); 115 penalty = @(v) -tau(v);
116 otherwise 116 otherwise
117 error('No such boundary condition: type = %s',type); 117 error('No such boundary condition: type = %s',type);
118 end 118 end
119 119
120 120
121 end 121 end
122 122
123 % Ruturns the boundary ops, half-norm, boundary indices and sign for the boundary specified by the string boundary. 123 % Returns the boundary sign. The right boundary is considered the positive boundary
124 % The right boundary for each coordinate direction is considered the positive boundary 124 % boundary -- string
125 function [e, H_b, index, s] = get_boundary_ops(obj, boundary) 125 function s = getBoundarySign(obj, boundary)
126 ind = grid.funcToMatrix(obj.grid, 1:obj.grid.N()); 126 assertIsMember(boundary, {'w', 'e', 's', 'n'})
127 switch boundary 127 switch boundary
128 case {'w', 'W', 'west', 'West'} 128 case {'e','n'}
129 e = obj.e_w; 129 s = 1;
130 H_b = obj.H_y; 130 case {'w','s'}
131 index = ind(1,:);
132 s = -1; 131 s = -1;
133 case {'e', 'E', 'east', 'East'} 132 end
134 e = obj.e_e; 133 end
135 H_b = obj.H_y; 134
136 index = ind(end,:); 135 % Returns the boundary operator op for the boundary specified by the string boundary.
137 s = 1; 136 % op -- string
138 case {'s', 'S', 'south', 'South'} 137 % boundary -- string
139 e = obj.e_s; 138 function o = getBoundaryOperator(obj, op, boundary)
140 H_b = obj.H_x; 139 assertIsMember(op, {'e'})
141 index = ind(:,1); 140 assertIsMember(boundary, {'w', 'e', 's', 'n'})
142 s = -1; 141
143 case {'n', 'N', 'north', 'North'} 142 o = obj.([op, '_', boundary]);
144 e = obj.e_n; 143 end
145 H_b = obj.H_x; 144
146 index = ind(:,end); 145 % Returns square boundary quadrature matrix, of dimension
147 s = 1; 146 % corresponding to the number of boundary points
148 otherwise 147 %
149 error('No such boundary: boundary = %s',boundary); 148 % boundary -- string
149 function H_b = getBoundaryQuadrature(obj, boundary)
150 assertIsMember(boundary, {'w', 'e', 's', 'n'})
151 H_b = obj.(['H_', boundary]);
152 end
153
154 % Returns square boundary quadrature matrix, of dimension
155 % corresponding to the number of boundary points
156 %
157 % boundary -- string
158 function H_1d = getOneDirectionalNorm(obj, boundary)
159 assertIsMember(boundary, {'w', 'e', 's', 'n'})
160 switch boundary
161 case {'w','e'}
162 H_1d = obj.H_y;
163 case {'s','n'}
164 H_1d = obj.H_x;
165 end
166 end
167
168 % Returns the indices of the boundary points in the grid matrix
169 % boundary -- string
170 function I = getBoundaryIndices(obj, boundary)
171 assertIsMember(boundary, {'w', 'e', 's', 'n'})
172
173 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.grid.m));
174 switch boundary
175 case 'w'
176 I = ind(1,:);
177 case 'e'
178 I = ind(end,:);
179 case 's'
180 I = ind(:,1)';
181 case 'n'
182 I = ind(:,end)';
150 end 183 end
151 end 184 end
152 185
153 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 186 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
154 error('An interface function does not exist yet'); 187 error('An interface function does not exist yet');