Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariableAnisotropic.m @ 1302:a0d615bde7f8 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add the hollow option to the anisotropic diffops
author Martin Almquist <malmquist@stanford.edu>
date Fri, 10 Jul 2020 20:24:23 -0700
parents 7fbee3de0ab2
children 49e3870335ef
comparison
equal deleted inserted replaced
1295:cb053fabbedc 1302:a0d615bde7f8
53 53
54 methods 54 methods
55 55
56 % The coefficients can either be function handles or grid functions 56 % The coefficients can either be function handles or grid functions
57 % optFlag -- if true, extra computations are performed, which may be helpful for optimization. 57 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
58 function obj = Elastic2dVariableAnisotropic(g, order, rho, C, opSet, optFlag) 58 function obj = Elastic2dVariableAnisotropic(g, order, rho, C, opSet, optFlag, hollow)
59 default_arg('hollow', false);
59 default_arg('rho', @(x,y) 0*x+1); 60 default_arg('rho', @(x,y) 0*x+1);
60 default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible}); 61 default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible});
61 default_arg('optFlag', false); 62 default_arg('optFlag', false);
62 dim = 2; 63 dim = 2;
63 64
194 195
195 % D2 196 % D2
196 switch order 197 switch order
197 case 2 198 case 2
198 width = 3; 199 width = 3;
200 nBP = 2;
199 case 4 201 case 4
200 width = 5; 202 width = 5;
203 nBP = 6;
201 case 6 204 case 6
202 width = 7; 205 width = 7;
206 nBP = 9;
203 end 207 end
204 for j = 1:dim 208 for j = 1:dim
205 for k = 1:dim 209 for k = 1:dim
206 for l = 1:dim 210 for l = 1:dim
207 D2_temp{j,k,l} = spalloc(m_tot, m_tot, width*m_tot); 211 if hollow
212 D2_temp{j,k,l} = sparse(m_tot, m_tot);
213 else
214 D2_temp{j,k,l} = spalloc(m_tot, m_tot, width*m_tot);
215 end
208 end 216 end
209 end 217 end
210 end 218 end
211 ind = grid.funcToMatrix(g, 1:m_tot); 219 ind = grid.funcToMatrix(g, 1:m_tot);
212 220
213 k = 1; 221 k = 1;
222 if hollow
223 mask = sparse(m(1), m(1));
224 mask(1:nBP, 1:nBP) = speye(nBP, nBP);
225 mask(end-nBP+1:end, end-nBP+1:end) = speye(nBP, nBP);
226 maskX = kron(mask, speye(m(2), m(2)));
227 maskX = E{1}*maskX*E{1}' + E{2}*maskX*E{2}';
228 end
214 for r = 1:m(2) 229 for r = 1:m(2)
215 p = ind(:,r); 230 p = ind(:,r);
216 for j = 1:dim 231 for j = 1:dim
217 for l = 1:dim 232 for l = 1:dim
218 coeff = C{k,j,k,l}; 233 coeff = C{k,j,k,l};
219 D_kk = D2{1}(coeff(p)); 234 D_kk = D2{1}(coeff(p));
235 if hollow && r > nBP && r < m(2) - nBP + 1
236 D_kk = mask*D_kk;
237 end
220 D2_temp{j,k,l}(p,p) = D_kk; 238 D2_temp{j,k,l}(p,p) = D_kk;
221 end 239 end
222 end 240 end
223 end 241 end
224 242
225 k = 2; 243 k = 2;
244 if hollow
245 mask = sparse(m(2), m(2));
246 mask(1:nBP, 1:nBP) = speye(nBP, nBP);
247 mask(end-nBP+1:end, end-nBP+1:end) = speye(nBP, nBP);
248 maskY = kron(speye(m(1), m(1)), mask);
249 maskY = E{1}*maskY*E{1}' + E{2}*maskY*E{2}';
250 end
226 for r = 1:m(1) 251 for r = 1:m(1)
227 p = ind(r,:); 252 p = ind(r,:);
228 for j = 1:dim 253 for j = 1:dim
229 for l = 1:dim 254 for l = 1:dim
230 coeff = C{k,j,k,l}; 255 coeff = C{k,j,k,l};
231 D_kk = D2{2}(coeff(p)); 256 D_kk = D2{2}(coeff(p));
257 if hollow && r > nBP && r < m(1) - nBP + 1
258 D_kk = mask*D_kk;
259 end
232 D2_temp{j,k,l}(p,p) = D_kk; 260 D2_temp{j,k,l}(p,p) = D_kk;
233 end 261 end
234 end 262 end
235 end 263 end
236 264
259 end 287 end
260 end 288 end
261 end 289 end
262 end 290 end
263 clear D2_temp; 291 clear D2_temp;
292 if hollow
293 mask = maskX + maskY;
294 mask = mask>0;
295
296 % D = maskX*maskY*D;
297 D = mask*D;
298 end
264 D = obj.RHOi_kron*D; 299 D = obj.RHOi_kron*D;
265 obj.D = D; 300 obj.D = D;
266 clear D; 301 clear D;
267 %=========================================%' 302 %=========================================%'
268 303