Mercurial > repos > public > sbplib

comparison +scheme/Heat2dVariable.m @ 905:459eeb99130f feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Include type as (optional) input parameter in the interface method of all schemes.
author Martin Almquist <malmquist@stanford.edu>
date Thu, 22 Nov 2018 22:03:44 -0800
parents 60eb7f46d8d9
children b9c98661ff5d
comparison
equal deleted inserted replaced
904:14b093a344eb 905:459eeb99130f
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
26 D2_kappa 26 D2_kappa
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 31
32 H_boundary % Boundary inner products 32 H_boundary % Boundary inner products
33 33
34 end 34 end
35 35
36 methods 36 methods
158 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 158 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
159 default_arg('type','Neumann'); 159 default_arg('type','Neumann');
160 default_arg('parameter', []); 160 default_arg('parameter', []);
161 161
162 % j is the coordinate direction of the boundary 162 % j is the coordinate direction of the boundary
163 % nj: outward unit normal component. 163 % nj: outward unit normal component.
164 % nj = -1 for west, south, bottom boundaries 164 % nj = -1 for west, south, bottom boundaries
165 % nj = 1 for east, north, top boundaries 165 % nj = 1 for east, north, top boundaries
166 [j, nj] = obj.get_boundary_number(boundary); 166 [j, nj] = obj.get_boundary_number(boundary);
167 switch nj 167 switch nj
168 case 1 168 case 1
174 end 174 end
175 175
176 Hi = obj.Hi; 176 Hi = obj.Hi;
177 H_gamma = obj.H_boundary{j}; 177 H_gamma = obj.H_boundary{j};
178 KAPPA = obj.KAPPA; 178 KAPPA = obj.KAPPA;
179 kappa_gamma = e{j}'*KAPPA*e{j}; 179 kappa_gamma = e{j}'*KAPPA*e{j};
180 180
181 switch type 181 switch type
182 182
183 % Dirichlet boundary condition 183 % Dirichlet boundary condition
184 case {'D','d','dirichlet','Dirichlet'} 184 case {'D','d','dirichlet','Dirichlet'}
185 closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' ); 185 closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' );
186 penalty = nj*Hi*d{j}*kappa_gamma*H_gamma; 186 penalty = nj*Hi*d{j}*kappa_gamma*H_gamma;
187 187
188 % Free boundary condition 188 % Free boundary condition
189 case {'N','n','neumann','Neumann'} 189 case {'N','n','neumann','Neumann'}
190 closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' ); 190 closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' );
191 penalty = nj*Hi*e{j}*kappa_gamma*H_gamma; 191 penalty = nj*Hi*e{j}*kappa_gamma*H_gamma;
192 192
193 % Unknown boundary condition 193 % Unknown boundary condition
194 otherwise 194 otherwise
195 error('No such boundary condition: type = %s',type); 195 error('No such boundary condition: type = %s',type);
196 end 196 end
197 end 197 end
198 198
199 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 199 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
200 % u denotes the solution in the own domain 200 % u denotes the solution in the own domain
201 % v denotes the solution in the neighbour domain 201 % v denotes the solution in the neighbour domain
202 error('Interface not implemented'); 202 error('Interface not implemented');
203 end 203 end
204 204