Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinearNewCorner.m @ 1117:27019aca2f17 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix bug in LaplCurveNewCorner that makes bc matrix non-sparse.
author Martin Almquist <malmquist@stanford.edu>
date Sat, 04 May 2019 14:44:34 -0700
parents 753de514ae77
children
comparison
equal deleted inserted replaced
1080:753de514ae77 1117:27019aca2f17
61 default_arg('b', @(x,y) 0*x + 1); 61 default_arg('b', @(x,y) 0*x + 1);
62 62
63 % assert(isa(g, 'grid.Curvilinear')) 63 % assert(isa(g, 'grid.Curvilinear'))
64 if isa(a, 'function_handle') 64 if isa(a, 'function_handle')
65 a = grid.evalOn(g, a); 65 a = grid.evalOn(g, a);
66 a = spdiag(a); 66 end
67 end 67 a = spdiag(a);
68 68
69 if isa(b, 'function_handle') 69 if isa(b, 'function_handle')
70 b = grid.evalOn(g, b); 70 b = grid.evalOn(g, b);
71 b = spdiag(b); 71 end
72 end 72 b = spdiag(b);
73 73
74 % If b is scalar 74 % If b is scalar
75 if length(b) == 1 75 if length(b) == 1
76 b = b*speye(g.N(), g.N()); 76 b = b*speye(g.N(), g.N());
77 end 77 end
303 switch type 303 switch type
304 % Dirichlet boundary condition 304 % Dirichlet boundary condition
305 case {'D','d','dirichlet'} 305 case {'D','d','dirichlet'}
306 tuning = 1.0; 306 tuning = 1.0;
307 307
308 sigma = 0; 308 sigma = 0*b_b;
309 for i = 1:obj.dim 309 for i = 1:obj.dim
310 sigma = sigma + e'*J*K{i,m}*K{i,m}*e; 310 sigma = sigma + e'*J*K{i,m}*K{i,m}*e;
311 end 311 end
312 sigma = sigma/s_b; 312 sigma = sigma/s_b;
313 % tau = tuning*(1/th_R + obj.dim/th_H)*sigma; 313 % tau = tuning*(1/th_R + obj.dim/th_H)*sigma;
399 K_v = v.K; 399 K_v = v.K;
400 J_v = v.J; 400 J_v = v.J;
401 b_b_v = e_v'*v.b*e_v; 401 b_b_v = e_v'*v.b*e_v;
402 402
403 %--- Penalty strength tau ------------- 403 %--- Penalty strength tau -------------
404 sigma_u = 0; 404 sigma_u = 0*b_b_u;
405 sigma_v = 0; 405 sigma_v = 0*b_b_v;
406 for i = 1:obj.dim 406 for i = 1:obj.dim
407 sigma_u = sigma_u + e_u'*J_u*K_u{i,m_u}*K_u{i,m_u}*e_u; 407 sigma_u = sigma_u + e_u'*J_u*K_u{i,m_u}*K_u{i,m_u}*e_u;
408 sigma_v = sigma_v + e_v'*J_v*K_v{i,m_v}*K_v{i,m_v}*e_v; 408 sigma_v = sigma_v + e_v'*J_v*K_v{i,m_v}*K_v{i,m_v}*e_v;
409 end 409 end
410 sigma_u = sigma_u/s_b_u; 410 sigma_u = sigma_u/s_b_u;