Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariableAnisotropic.m @ 1253:89dad61cad22 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make Elastic2dVariable faster and more memory efficient
author Martin Almquist <malmquist@stanford.edu>
date Tue, 04 Feb 2020 10:15:42 -0800
parents 8fc2f9a4c882
children 7fbee3de0ab2
comparison
equal deleted inserted replaced
1252:8fc2f9a4c882 1253:89dad61cad22
19 RHO, RHOi, RHOi_kron % Density 19 RHO, RHOi, RHOi_kron % Density
20 C % Elastic stiffness tensor 20 C % Elastic stiffness tensor
21 21
22 D % Total operator 22 D % Total operator
23 D1 % First derivatives 23 D1 % First derivatives
24 D2 % Second derivatives 24 % D2 % Second derivatives
25 25
26 % Boundary operators in cell format, used for BC 26 % Boundary operators in cell format, used for BC
27 T_w, T_e, T_s, T_n 27 T_w, T_e, T_s, T_n
28 28
29 % Traction operators 29 % Traction operators
145 obj.RHO = RHO; 145 obj.RHO = RHO;
146 obj.RHOi = inv(RHO); 146 obj.RHOi = inv(RHO);
147 obj.RHOi_kron = kron(obj.RHOi, I_dim); 147 obj.RHOi_kron = kron(obj.RHOi, I_dim);
148 148
149 obj.D1 = cell(dim,1); 149 obj.D1 = cell(dim,1);
150 obj.D2 = cell(dim,dim,dim); 150 D2_temp = cell(dim,dim,dim);
151 151
152 % D1 152 % D1
153 obj.D1{1} = kron(D1{1},I{2}); 153 obj.D1{1} = kron(D1{1},I{2});
154 obj.D1{2} = kron(I{1},D1{2}); 154 obj.D1{2} = kron(I{1},D1{2});
155 155
191 e2_s = (e_scalar_s'*E{2}')'; 191 e2_s = (e_scalar_s'*E{2}')';
192 e2_n = (e_scalar_n'*E{2}')'; 192 e2_n = (e_scalar_n'*E{2}')';
193 193
194 194
195 % D2 195 % D2
196 switch order
197 case 2
198 width = 3;
199 case 4
200 width = 5;
201 case 6
202 width = 7;
203 end
196 for j = 1:dim 204 for j = 1:dim
197 for k = 1:dim 205 for k = 1:dim
198 for l = 1:dim 206 for l = 1:dim
199 obj.D2{j,k,l} = sparse(m_tot,m_tot); 207 D2_temp{j,k,l} = spalloc(m_tot, m_tot, width*m_tot);
200 end 208 end
201 end 209 end
202 end 210 end
203 ind = grid.funcToMatrix(g, 1:m_tot); 211 ind = grid.funcToMatrix(g, 1:m_tot);
204 212
207 p = ind(:,r); 215 p = ind(:,r);
208 for j = 1:dim 216 for j = 1:dim
209 for l = 1:dim 217 for l = 1:dim
210 coeff = C{k,j,k,l}; 218 coeff = C{k,j,k,l};
211 D_kk = D2{1}(coeff(p)); 219 D_kk = D2{1}(coeff(p));
212 obj.D2{j,k,l}(p,p) = D_kk; 220 D2_temp{j,k,l}(p,p) = D_kk;
213 end 221 end
214 end 222 end
215 end 223 end
216 224
217 k = 2; 225 k = 2;
219 p = ind(r,:); 227 p = ind(r,:);
220 for j = 1:dim 228 for j = 1:dim
221 for l = 1:dim 229 for l = 1:dim
222 coeff = C{k,j,k,l}; 230 coeff = C{k,j,k,l};
223 D_kk = D2{2}(coeff(p)); 231 D_kk = D2{2}(coeff(p));
224 obj.D2{j,k,l}(p,p) = D_kk; 232 D2_temp{j,k,l}(p,p) = D_kk;
225 end 233 end
226 end 234 end
227 end 235 end
228 236
229 % Quadratures 237 % Quadratures
234 obj.H_s = H{1}; 242 obj.H_s = H{1};
235 obj.H_n = H{1}; 243 obj.H_n = H{1};
236 obj.H_1D = {H{1}, H{2}}; 244 obj.H_1D = {H{1}, H{2}};
237 245
238 % Differentiation matrix D (without SAT) 246 % Differentiation matrix D (without SAT)
239 D2 = obj.D2; 247 D2_temp;
240 D1 = obj.D1; 248 D1 = obj.D1;
241 D = sparse(dim*m_tot,dim*m_tot); 249 D = sparse(dim*m_tot,dim*m_tot);
242 for i = 1:dim 250 for i = 1:dim
243 for j = 1:dim 251 for j = 1:dim
244 for k = 1:dim 252 for k = 1:dim
245 for l = 1:dim 253 for l = 1:dim
246 if i == k 254 if i == k
247 D = D + E{j}*D2{j,k,l}*E{l}'; 255 D = D + E{j}*D2_temp{j,k,l}*E{l}';
256 D2_temp{j,k,l} = [];
248 else 257 else
249 D = D + E{j}*D1{i}*C_mat{i,j,k,l}*D1{k}*E{l}'; 258 D = D + E{j}*D1{i}*C_mat{i,j,k,l}*D1{k}*E{l}';
250 end 259 end
251 end 260 end
252 end 261 end
253 end 262 end
254 end 263 end
264 clear D2_temp;
255 D = obj.RHOi_kron*D; 265 D = obj.RHOi_kron*D;
256 obj.D = D; 266 obj.D = D;
267 clear D;
257 %=========================================%' 268 %=========================================%'
258 269
259 % Numerical traction operators for BC. 270 % Numerical traction operators for BC.
260 % 271 %
261 % Formula at boundary j: % tau^{j}_i = sum_l T^{j}_{il} u_l 272 % Formula at boundary j: % tau^{j}_i = sum_l T^{j}_{il} u_l