Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariableAnisotropic.m @ 1206:6b203030fb37 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix performance (full matrices) issues in ElasticVariableAnisotropic
author Martin Almquist <malmquist@stanford.edu>
date Sun, 15 Sep 2019 14:39:42 -0700
parents 687515778437
children 7f427275bc9c
comparison
equal deleted inserted replaced
1205:3258dca12af8 1206:6b203030fb37
138 d1_0{i} = ops{i}.d1_l; 138 d1_0{i} = ops{i}.d1_l;
139 d1_m{i} = ops{i}.d1_r; 139 d1_m{i} = ops{i}.d1_r;
140 end 140 end
141 141
142 %====== Assemble full operators ======== 142 %====== Assemble full operators ========
143 I_dim = speye(dim, dim);
143 RHO = spdiag(rho); 144 RHO = spdiag(rho);
144 obj.RHO = RHO; 145 obj.RHO = RHO;
145 obj.RHOi = inv(RHO); 146 obj.RHOi = inv(RHO);
147 obj.RHOi_kron = kron(obj.RHOi, I_dim);
146 148
147 obj.D1 = cell(dim,1); 149 obj.D1 = cell(dim,1);
148 obj.D2 = cell(dim,dim,dim); 150 obj.D2 = cell(dim,dim,dim);
149 151
150 % D1 152 % D1
162 e_scalar_w = e_l{1}; 164 e_scalar_w = e_l{1};
163 e_scalar_e = e_r{1}; 165 e_scalar_e = e_r{1};
164 e_scalar_s = e_l{2}; 166 e_scalar_s = e_l{2};
165 e_scalar_n = e_r{2}; 167 e_scalar_n = e_r{2};
166 168
167 I_dim = speye(dim, dim);
168 e_w = kron(e_scalar_w, I_dim); 169 e_w = kron(e_scalar_w, I_dim);
169 e_e = kron(e_scalar_e, I_dim); 170 e_e = kron(e_scalar_e, I_dim);
170 e_s = kron(e_scalar_s, I_dim); 171 e_s = kron(e_scalar_s, I_dim);
171 e_n = kron(e_scalar_n, I_dim); 172 e_n = kron(e_scalar_n, I_dim);
172 173
191 e2_n = (e_scalar_n'*E{2}')'; 192 e2_n = (e_scalar_n'*E{2}')';
192 193
193 194
194 % D2 195 % D2
195 for i = 1:dim 196 for i = 1:dim
196 for k = 2:dim 197 for k = 1:dim
197 for l = 2:dim 198 for l = 1:dim
198 obj.D2{i,k,l} = sparse(m_tot); 199 obj.D2{i,k,l} = sparse(m_tot,m_tot);
199 end 200 end
200 end 201 end
201 end 202 end
202 ind = grid.funcToMatrix(g, 1:m_tot); 203 ind = grid.funcToMatrix(g, 1:m_tot);
203 204
236 237
237 % Differentiation matrix D (without SAT) 238 % Differentiation matrix D (without SAT)
238 D2 = obj.D2; 239 D2 = obj.D2;
239 D1 = obj.D1; 240 D1 = obj.D1;
240 D = sparse(dim*m_tot,dim*m_tot); 241 D = sparse(dim*m_tot,dim*m_tot);
241 d = @kroneckerDelta; % Kronecker delta
242 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
243 for i = 1:dim 242 for i = 1:dim
244 for j = 1:dim 243 for j = 1:dim
245 for k = 1:dim 244 for k = 1:dim
246 for l = 1:dim 245 for l = 1:dim
247 D = D + E{i}*inv(RHO)*( d(j,k)*D2{i,k,l}*E{l}' +... 246 if j == k
248 db(j,k)*D1{j}*C_mat{i,j,k,l}*D1{k}*E{l}' ... 247 D = D + E{i}*D2{i,k,l}*E{l}';
249 ); 248 else
249 D = D + E{i}*D1{j}*C_mat{i,j,k,l}*D1{k}*E{l}';
250 end
250 end 251 end
251 end 252 end
252 end 253 end
253 end 254 end
255 D = obj.RHOi_kron*D;
254 obj.D = D; 256 obj.D = D;
255 %=========================================%' 257 %=========================================%'
256 258
257 % Numerical traction operators for BC. 259 % Numerical traction operators for BC.
258 % 260 %
342 obj.tau_e = (e_e'*e1_e*obj.tau1_e')' + (e_e'*e2_e*obj.tau2_e')'; 344 obj.tau_e = (e_e'*e1_e*obj.tau1_e')' + (e_e'*e2_e*obj.tau2_e')';
343 obj.tau_s = (e_s'*e1_s*obj.tau1_s')' + (e_s'*e2_s*obj.tau2_s')'; 345 obj.tau_s = (e_s'*e1_s*obj.tau1_s')' + (e_s'*e2_s*obj.tau2_s')';
344 obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')'; 346 obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')';
345 347
346 % Kroneckered norms and coefficients 348 % Kroneckered norms and coefficients
347 obj.RHOi_kron = kron(obj.RHOi, I_dim);
348 obj.Hi_kron = kron(obj.Hi, I_dim); 349 obj.Hi_kron = kron(obj.Hi, I_dim);
349 350
350 % Misc. 351 % Misc.
351 obj.m = m; 352 obj.m = m;
352 obj.h = h; 353 obj.h = h;