Mercurial > repos > public > sbplib

comparison +scheme/Heat2dVariable.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 8d73fcdb07a5
children
comparison
equal deleted inserted replaced
1196:f6c571d8f22f 1197:433c89bf19e0
1 classdef Heat2dVariable < scheme.Scheme 1 classdef Heat2dVariable < scheme.Scheme
2 2
3 % Discretizes the Laplacian with variable coefficent, 3 % Discretizes the Laplacian with variable coefficent,
4 % In the Heat equation way (i.e., the discretization matrix is not necessarily 4 % In the Heat equation way (i.e., the discretization matrix is not necessarily
5 % symmetric) 5 % symmetric)
6 % u_t = div * (kappa * grad u ) 6 % u_t = div * (kappa * grad u )
7 % opSet should be cell array of opSets, one per dimension. This 7 % opSet should be cell array of opSets, one per dimension. This
8 % is useful if we have periodic BC in one direction. 8 % is useful if we have periodic BC in one direction.
9 9
10 properties 10 properties
11 m % Number of points in each direction, possibly a vector 11 m % Number of points in each direction, possibly a vector
27 27
28 H, Hi % Inner products 28 H, Hi % Inner products
29 e_l, e_r 29 e_l, e_r
30 d1_l, d1_r % Normal derivatives at the boundary 30 d1_l, d1_r % Normal derivatives at the boundary
31 alpha % Vector of borrowing constants 31 alpha % Vector of borrowing constants
32 32
33 H_boundary % Boundary inner products 33 H_boundary % Boundary inner products
34 34
35 end 35 end
36 36
37 methods 37 methods
160 function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning) 160 function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning)
161 default_arg('type','Neumann'); 161 default_arg('type','Neumann');
162 default_arg('symmetric', false); 162 default_arg('symmetric', false);
163 default_arg('tuning',1.2); 163 default_arg('tuning',1.2);
164 164
165 % j is the coordinate direction of the boundary 165 % nj: outward unit normal component.
166 % nj: outward unit normal component.
167 % nj = -1 for west, south, bottom boundaries 166 % nj = -1 for west, south, bottom boundaries
168 % nj = 1 for east, north, top boundaries 167 % nj = 1 for east, north, top boundaries
169 [j, nj] = obj.get_boundary_number(boundary); 168 nj = obj.getBoundarySign(boundary);
170 switch nj
171 case 1
172 e = obj.e_r;
173 d = obj.d1_r;
174 case -1
175 e = obj.e_l;
176 d = obj.d1_l;
177 end
178 169
179 Hi = obj.Hi; 170 Hi = obj.Hi;
180 H_gamma = obj.H_boundary{j}; 171 [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary);
172 H_gamma = obj.getBoundaryQuadrature(boundary);
173 alpha = obj.getBoundaryBorrowing(boundary);
174
181 KAPPA = obj.KAPPA; 175 KAPPA = obj.KAPPA;
182 kappa_gamma = e{j}'*KAPPA*e{j}; 176 kappa_gamma = e'*KAPPA*e;
183 h = obj.h(j);
184 alpha = h*obj.alpha(j);
185 177
186 switch type 178 switch type
187 179
188 % Dirichlet boundary condition 180 % Dirichlet boundary condition
189 case {'D','d','dirichlet','Dirichlet'} 181 case {'D','d','dirichlet','Dirichlet'}
190 182
191 if ~symmetric 183 if ~symmetric
192 closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' ); 184 closure = -nj*Hi*d*kappa_gamma*H_gamma*(e' );
193 penalty = nj*Hi*d{j}*kappa_gamma*H_gamma; 185 penalty = nj*Hi*d*kappa_gamma*H_gamma;
194 else 186 else
195 closure = nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' )... 187 closure = nj*Hi*d*kappa_gamma*H_gamma*(e' )...
196 -tuning*2/alpha*Hi*e{j}*kappa_gamma*H_gamma*(e{j}' ) ; 188 -tuning*2/alpha*Hi*e*kappa_gamma*H_gamma*(e' ) ;
197 penalty = -nj*Hi*d{j}*kappa_gamma*H_gamma ... 189 penalty = -nj*Hi*d*kappa_gamma*H_gamma ...
198 +tuning*2/alpha*Hi*e{j}*kappa_gamma*H_gamma; 190 +tuning*2/alpha*Hi*e*kappa_gamma*H_gamma;
199 end 191 end
200 192
201 % Free boundary condition 193 % Free boundary condition
202 case {'N','n','neumann','Neumann'} 194 case {'N','n','neumann','Neumann'}
203 closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' ); 195 closure = -nj*Hi*e*kappa_gamma*H_gamma*(d' );
204 penalty = Hi*e{j}*kappa_gamma*H_gamma; 196 penalty = Hi*e*kappa_gamma*H_gamma;
205 % penalty is for normal derivative and not for derivative, hence the sign. 197 % penalty is for normal derivative and not for derivative, hence the sign.
206 198
207 % Unknown boundary condition 199 % Unknown boundary condition
208 otherwise 200 otherwise
209 error('No such boundary condition: type = %s',type); 201 error('No such boundary condition: type = %s',type);
214 % u denotes the solution in the own domain 206 % u denotes the solution in the own domain
215 % v denotes the solution in the neighbour domain 207 % v denotes the solution in the neighbour domain
216 error('Interface not implemented'); 208 error('Interface not implemented');
217 end 209 end
218 210
219 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 211 % Returns the boundary operator op for the boundary specified by the string boundary.
220 function [j, nj] = get_boundary_number(obj, boundary) 212 % op -- string or a cell array of strings
221 213 % boundary -- string
222 switch boundary 214 function varargout = getBoundaryOperator(obj, op, boundary)
223 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 215 assertIsMember(boundary, {'w', 'e', 's', 'n'})
224 j = 1; 216
225 case {'s','S','south','South', 'n', 'N', 'north', 'North'} 217 if ~iscell(op)
226 j = 2; 218 op = {op};
227 otherwise 219 end
228 error('No such boundary: boundary = %s',boundary); 220
229 end 221 for i = 1:numel(op)
230 222 switch op{i}
231 switch boundary
232 case {'w','W','west','West','s','S','south','South'}
233 nj = -1;
234 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
235 nj = 1;
236 end
237 end
238
239 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
240 function [return_op] = get_boundary_operator(obj, op, boundary)
241
242 switch boundary
243 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
244 j = 1;
245 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
246 j = 2;
247 otherwise
248 error('No such boundary: boundary = %s',boundary);
249 end
250
251 switch op
252 case 'e' 223 case 'e'
253 switch boundary 224 switch boundary
254 case {'w','W','west','West','s','S','south','South'} 225 case 'w'
255 return_op = obj.e_l{j}; 226 e = obj.e_l{1};
256 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 227 case 'e'
257 return_op = obj.e_r{j}; 228 e = obj.e_r{1};
229 case 's'
230 e = obj.e_l{2};
231 case 'n'
232 e = obj.e_r{2};
258 end 233 end
234 varargout{i} = e;
235
259 case 'd' 236 case 'd'
260 switch boundary 237 switch boundary
261 case {'w','W','west','West','s','S','south','South'} 238 case 'w'
262 return_op = obj.d1_l{j}; 239 d = obj.d1_l{1};
263 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 240 case 'e'
264 return_op = obj.d1_r{j}; 241 d = obj.d1_r{1};
242 case 's'
243 d = obj.d1_l{2};
244 case 'n'
245 d = obj.d1_r{2};
265 end 246 end
266 otherwise 247 varargout{i} = d;
267 error(['No such operator: operatr = ' op]); 248 end
268 end 249 end
269 250 end
251
252 % Returns square boundary quadrature matrix, of dimension
253 % corresponding to the number of boundary points
254 %
255 % boundary -- string
256 function H_b = getBoundaryQuadrature(obj, boundary)
257 assertIsMember(boundary, {'w', 'e', 's', 'n'})
258
259 switch boundary
260 case 'w'
261 H_b = obj.H_boundary{1};
262 case 'e'
263 H_b = obj.H_boundary{1};
264 case 's'
265 H_b = obj.H_boundary{2};
266 case 'n'
267 H_b = obj.H_boundary{2};
268 end
269 end
270
271 % Returns the boundary sign. The right boundary is considered the positive boundary
272 % boundary -- string
273 function s = getBoundarySign(obj, boundary)
274 assertIsMember(boundary, {'w', 'e', 's', 'n'})
275
276 switch boundary
277 case {'e','n'}
278 s = 1;
279 case {'w','s'}
280 s = -1;
281 end
282 end
283
284 % Returns borrowing constant gamma*h
285 % boundary -- string
286 function gamm = getBoundaryBorrowing(obj, boundary)
287 assertIsMember(boundary, {'w', 'e', 's', 'n'})
288
289 switch boundary
290 case {'w','e'}
291 gamm = obj.h(1)*obj.alpha(1);
292 case {'s','n'}
293 gamm = obj.h(2)*obj.alpha(2);
294 end
270 end 295 end
271 296
272 function N = size(obj) 297 function N = size(obj)
273 N = prod(obj.m); 298 N = prod(obj.m);
274 end 299 end