Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 551:c5a7a13c03dc feature/grids/laplace_refactor

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Switch to using spdiag instead of spdiags
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 28 Aug 2017 19:01:40 +0200
parents 08b6281ba2a9
children ffaf13533c27
comparison
equal deleted inserted replaced
550:e860670e72f1 551:c5a7a13c03dc
103 a12 = -1./J .* (x_u.*x_v + y_u.*y_v); 103 a12 = -1./J .* (x_u.*x_v + y_u.*y_v);
104 a22 = 1./J .* (x_u.^2 + y_u.^2); 104 a22 = 1./J .* (x_u.^2 + y_u.^2);
105 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2)); 105 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2));
106 106
107 % Assemble full operators 107 % Assemble full operators
108 L_12 = spdiags(a12, 0, m_tot, m_tot); 108 L_12 = spdiag(a12);
109 Duv = Du*L_12*Dv; 109 Duv = Du*L_12*Dv;
110 Dvu = Dv*L_12*Du; 110 Dvu = Dv*L_12*Du;
111 111
112 Duu = sparse(m_tot); 112 Duu = sparse(m_tot);
113 Dvv = sparse(m_tot); 113 Dvv = sparse(m_tot);
155 obj.order = order; 155 obj.order = order;
156 obj.grid = g; 156 obj.grid = g;
157 157
158 obj.a = a; 158 obj.a = a;
159 obj.b = b; 159 obj.b = b;
160 obj.J = spdiags(J, 0, m_tot, m_tot); 160 obj.J = spdiag(J);
161 obj.Ji = spdiags(1./J, 0, m_tot, m_tot); 161 obj.Ji = spdiag(1./J);
162 obj.a11 = a11; 162 obj.a11 = a11;
163 obj.a12 = a12; 163 obj.a12 = a12;
164 obj.a22 = a22; 164 obj.a22 = a22;
165 obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv); 165 obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv);
166 obj.lambda = lambda; 166 obj.lambda = lambda;
213 penalty = -obj.Ji*obj.a * penalty_parameter_1; 213 penalty = -obj.Ji*obj.a * penalty_parameter_1;
214 214
215 215
216 % Neumann boundary condition 216 % Neumann boundary condition
217 case {'N','n','neumann'} 217 case {'N','n','neumann'}
218 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); 218 a_n = spdiag(coeff_n);
219 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); 219 a_t = spdiag(coeff_t);
220 d = (a_n * d_n' + a_t*d_t')'; 220 d = (a_n * d_n' + a_t*d_t')';
221 221
222 tau1 = -s; 222 tau1 = -s;
223 tau2 = 0; 223 tau2 = 0;
224 tau = obj.a * obj.Ji*(tau1*e + tau2*d); 224 tau = obj.a * obj.Ji*(tau1*e + tau2*d);
229 % Characteristic boundary condition 229 % Characteristic boundary condition
230 case {'characteristic', 'char', 'c'} 230 case {'characteristic', 'char', 'c'}
231 default_arg('parameter', 1); 231 default_arg('parameter', 1);
232 beta = parameter; 232 beta = parameter;
233 233
234 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); 234 a_n = spdiag(coeff_n);
235 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); 235 a_t = spdiag(coeff_t);
236 d = s*(a_n * d_n' + a_t*d_t')'; % outward facing normal derivative 236 d = s*(a_n * d_n' + a_t*d_t')'; % outward facing normal derivative
237 237
238 tau = -obj.a * 1/beta*obj.Ji*e; 238 tau = -obj.a * 1/beta*obj.Ji*e;
239 239
240 closure{1} = halfnorm_inv*tau*spdiag(scale_factor)*e'; 240 closure{1} = halfnorm_inv*tau*spdiag(scale_factor)*e';