Mercurial > repos > public > sbplib

comparison +scheme/Heat2dVariable.m @ 725:e9e15d64f9f9 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement symmetric Dirichlet BC with variable tuning parameter in Heat2D scheme.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 30 Mar 2018 15:35:46 -0700
parents a9e8c9d71307
children 21394c78c72e
comparison
equal deleted inserted replaced
724:a9e8c9d71307 725:e9e15d64f9f9
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 alpha % Vector of borrowing constants
31 32
32 H_boundary % Boundary inner products 33 H_boundary % Boundary inner products
33 34
34 end 35 end
35 36
142 obj.m = m; 143 obj.m = m;
143 obj.h = h; 144 obj.h = h;
144 obj.order = order; 145 obj.order = order;
145 obj.grid = g; 146 obj.grid = g;
146 obj.dim = dim; 147 obj.dim = dim;
148 obj.alpha = [ops{1}.borrowing.M.d1, ops{2}.borrowing.M.d1];
147 149
148 end 150 end
149 151
150 152
151 % Closure functions return the operators applied to the own domain to close the boundary 153 % Closure functions return the operators applied to the own domain to close the boundary
153 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 155 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
154 % type is a string specifying the type of boundary condition. 156 % type is a string specifying the type of boundary condition.
155 % data is a function returning the data that should be applied at the boundary. 157 % data is a function returning the data that should be applied at the boundary.
156 % neighbour_scheme is an instance of Scheme that should be interfaced to. 158 % neighbour_scheme is an instance of Scheme that should be interfaced to.
157 % neighbour_boundary is a string specifying which boundary to interface to. 159 % neighbour_boundary is a string specifying which boundary to interface to.
158 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 160 function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning)
159 default_arg('type','Neumann'); 161 default_arg('type','Neumann');
160 default_arg('parameter', []); 162 default_arg('symmetric', false);
163 default_arg('tuning',1.2);
161 164
162 % j is the coordinate direction of the boundary 165 % j is the coordinate direction of the boundary
163 % nj: outward unit normal component. 166 % nj: outward unit normal component.
164 % nj = -1 for west, south, bottom boundaries 167 % nj = -1 for west, south, bottom boundaries
165 % nj = 1 for east, north, top boundaries 168 % nj = 1 for east, north, top boundaries
175 178
176 Hi = obj.Hi; 179 Hi = obj.Hi;
177 H_gamma = obj.H_boundary{j}; 180 H_gamma = obj.H_boundary{j};
178 KAPPA = obj.KAPPA; 181 KAPPA = obj.KAPPA;
179 kappa_gamma = e{j}'*KAPPA*e{j}; 182 kappa_gamma = e{j}'*KAPPA*e{j};
183 h = obj.h(j);
184 alpha = h*obj.alpha(j);
180 185
181 switch type 186 switch type
182 187
183 % Dirichlet boundary condition 188 % Dirichlet boundary condition
184 case {'D','d','dirichlet','Dirichlet'} 189 case {'D','d','dirichlet','Dirichlet'}
190
191 if ~symmetric
185 closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' ); 192 closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' );
186 penalty = nj*Hi*d{j}*kappa_gamma*H_gamma; 193 penalty = nj*Hi*d{j}*kappa_gamma*H_gamma;
194 else
195 closure = nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' )...
196 -tuning*2/alpha*Hi*e{j}*kappa_gamma*H_gamma*(e{j}' ) ;
197 penalty = -nj*Hi*d{j}*kappa_gamma*H_gamma ...
198 +tuning*2/alpha*Hi*e{j}*kappa_gamma*H_gamma;
199 end
187 200
188 % Free boundary condition 201 % Free boundary condition
189 case {'N','n','neumann','Neumann'} 202 case {'N','n','neumann','Neumann'}
190 closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' ); 203 closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' );
191 penalty = Hi*e{j}*kappa_gamma*H_gamma; 204 penalty = Hi*e{j}*kappa_gamma*H_gamma;