Mercurial > repos > public > sbplib

comparison +scheme/elasticShearVariable.m @ 667:ed853945ee99 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make free BC work with data. Bugfix domain size.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 20 Dec 2017 06:52:17 +0100
parents 8e6dfd22fc59
children 17e62551cdc2
comparison
equal deleted inserted replaced
664:8e6dfd22fc59 667:ed853945ee99
42 rho = grid.evalOn(g, rho_fun); 42 rho = grid.evalOn(g, rho_fun);
43 m = g.size(); 43 m = g.size();
44 m_tot = g.N(); 44 m_tot = g.N();
45 45
46 h = g.scaling(); 46 h = g.scaling();
47 L = (m-1).*h;
47 48
48 % 1D operators 49 % 1D operators
49 ops = cell(dim,1); 50 ops = cell(dim,1);
50 for i = 1:dim 51 for i = 1:dim
51 ops{i} = opSet(m(i), {0, 1}, order); 52 ops{i} = opSet(m(i), {0, L(i)}, order);
52 end 53 end
53 54
54 I = cell(dim,1); 55 I = cell(dim,1);
55 D1 = cell(dim,1); 56 D1 = cell(dim,1);
56 D2 = cell(dim,1); 57 D2 = cell(dim,1);
228 penalty = -obj.a*obj.Hi*tau; 229 penalty = -obj.a*obj.Hi*tau;
229 230
230 231
231 % Free boundary condition 232 % Free boundary condition
232 case {'F','f','Free','free'} 233 case {'F','f','Free','free'}
233 closure = sparse(obj.dim*obj.grid.N,obj.dim*obj.grid.N); 234 closures = cell(obj.dim,1);
234 penalty = sparse(obj.dim*obj.grid.N,obj.dim*obj.grid.N); 235 penalties = cell(obj.dim,1);
235 % Loop over components 236 % Loop over components
236 for i = 1:obj.dim 237 for i = 1:obj.dim
237 closure = closure + E{i}*inv(RHO)*(-nj)*Hi*e{j}*H_gamma*(... 238 closures{i} = E{i}*inv(RHO)*(-nj)*Hi*e{j}*H_gamma*(...
238 e{j}'*A*e{j}*d{j}'*E{i}' + ... 239 e{j}'*A*e{j}*d{j}'*E{i}' + ...
239 delta(i,j)*e{j}'*A*e{i}*d{i}'*E{j}' + ... 240 delta(i,j)*e{j}'*A*e{i}*d{i}'*E{j}' + ...
240 delta_b(i,j)*e{j}'*A*D1{i}*E{j}' ... 241 delta_b(i,j)*e{j}'*A*D1{i}*E{j}' ...
241 ); 242 );
242 penalty = penalty - E{i}*inv(RHO)*(-nj)*Hi*e{j}*H_gamma*e{j}'*E{j}'; 243 penalties{i} = - E{i}*inv(RHO)*(-nj)*Hi*e{j}*H_gamma;
243 end 244 end
245 [rows, cols] = size(closures{1});
246 closure = sparse(rows, cols);
247 for i = 1:obj.dim
248 closure = closure + closures{i};
249 end
250 penalty = penalties;
244 251
245 252
246 % Unknown boundary condition 253 % Unknown boundary condition
247 otherwise 254 otherwise
248 error('No such boundary condition: type = %s',type); 255 error('No such boundary condition: type = %s',type);