Mercurial > repos > public > sbplib
changeset 1302:a0d615bde7f8 feature/poroelastic
Add the hollow option to the anisotropic diffops
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Fri, 10 Jul 2020 20:24:23 -0700 |
| parents | cb053fabbedc |
| children | 49e3870335ef |
| files | +scheme/Elastic2dCurvilinearAnisotropic.m +scheme/Elastic2dVariableAnisotropic.m |
| diffstat | 2 files changed, 40 insertions(+), 4 deletions(-) [+] |
line wrap: on
line diff
diff -r cb053fabbedc -r a0d615bde7f8 +scheme/Elastic2dCurvilinearAnisotropic.m --- a/+scheme/Elastic2dCurvilinearAnisotropic.m Thu Jul 02 20:05:50 2020 -0700 +++ b/+scheme/Elastic2dCurvilinearAnisotropic.m Fri Jul 10 20:24:23 2020 -0700 @@ -66,7 +66,8 @@ % The coefficients can either be function handles or grid functions % optFlag -- if true, extra computations are performed, which may be helpful for optimization. - function obj = Elastic2dCurvilinearAnisotropic(g, order, rho, C, opSet, optFlag) + function obj = Elastic2dCurvilinearAnisotropic(g, order, rho, C, opSet, optFlag, hollow) + default_arg('hollow', false); default_arg('rho', @(x,y) 0*x+1); default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible}); default_arg('optFlag', false); @@ -188,7 +189,7 @@ end gRef = grid.equidistant([m_u, m_v], {0,1}, {0,1}); - refObj = scheme.Elastic2dVariableAnisotropic(gRef, order, rho_tilde, PHI, opSet); + refObj = scheme.Elastic2dVariableAnisotropic(gRef, order, rho_tilde, PHI, opSet, [], hollow); %---- Set object properties ------ obj.RHO = spdiag(rho);
diff -r cb053fabbedc -r a0d615bde7f8 +scheme/Elastic2dVariableAnisotropic.m --- a/+scheme/Elastic2dVariableAnisotropic.m Thu Jul 02 20:05:50 2020 -0700 +++ b/+scheme/Elastic2dVariableAnisotropic.m Fri Jul 10 20:24:23 2020 -0700 @@ -55,7 +55,8 @@ % The coefficients can either be function handles or grid functions % optFlag -- if true, extra computations are performed, which may be helpful for optimization. - function obj = Elastic2dVariableAnisotropic(g, order, rho, C, opSet, optFlag) + function obj = Elastic2dVariableAnisotropic(g, order, rho, C, opSet, optFlag, hollow) + default_arg('hollow', false); default_arg('rho', @(x,y) 0*x+1); default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible}); default_arg('optFlag', false); @@ -196,39 +197,66 @@ switch order case 2 width = 3; + nBP = 2; case 4 width = 5; + nBP = 6; case 6 width = 7; + nBP = 9; end for j = 1:dim for k = 1:dim for l = 1:dim - D2_temp{j,k,l} = spalloc(m_tot, m_tot, width*m_tot); + if hollow + D2_temp{j,k,l} = sparse(m_tot, m_tot); + else + D2_temp{j,k,l} = spalloc(m_tot, m_tot, width*m_tot); + end end end end ind = grid.funcToMatrix(g, 1:m_tot); k = 1; + if hollow + mask = sparse(m(1), m(1)); + mask(1:nBP, 1:nBP) = speye(nBP, nBP); + mask(end-nBP+1:end, end-nBP+1:end) = speye(nBP, nBP); + maskX = kron(mask, speye(m(2), m(2))); + maskX = E{1}*maskX*E{1}' + E{2}*maskX*E{2}'; + end for r = 1:m(2) p = ind(:,r); for j = 1:dim for l = 1:dim coeff = C{k,j,k,l}; D_kk = D2{1}(coeff(p)); + if hollow && r > nBP && r < m(2) - nBP + 1 + D_kk = mask*D_kk; + end D2_temp{j,k,l}(p,p) = D_kk; end end end k = 2; + if hollow + mask = sparse(m(2), m(2)); + mask(1:nBP, 1:nBP) = speye(nBP, nBP); + mask(end-nBP+1:end, end-nBP+1:end) = speye(nBP, nBP); + maskY = kron(speye(m(1), m(1)), mask); + maskY = E{1}*maskY*E{1}' + E{2}*maskY*E{2}'; + end for r = 1:m(1) p = ind(r,:); for j = 1:dim for l = 1:dim coeff = C{k,j,k,l}; D_kk = D2{2}(coeff(p)); + if hollow && r > nBP && r < m(1) - nBP + 1 + D_kk = mask*D_kk; + end D2_temp{j,k,l}(p,p) = D_kk; end end @@ -261,6 +289,13 @@ end end clear D2_temp; + if hollow + mask = maskX + maskY; + mask = mask>0; + + % D = maskX*maskY*D; + D = mask*D; + end D = obj.RHOi_kron*D; obj.D = D; clear D;