diff +scheme/Elastic2dVariableAnisotropic.m @ 1304:a38e80fdbf60 feature/poroelastic

Improve performance for the hollow case
author Martin Almquist <malmquist@stanford.edu>
date Sat, 11 Jul 2020 07:25:06 -0700
parents 49e3870335ef
children 8d9fc7981796
line wrap: on
line diff
--- a/+scheme/Elastic2dVariableAnisotropic.m	Sat Jul 11 06:54:15 2020 -0700
+++ b/+scheme/Elastic2dVariableAnisotropic.m	Sat Jul 11 07:25:06 2020 -0700
@@ -231,8 +231,8 @@
                 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}';
+                maskXSmall = kron(mask, speye(m(2), m(2)));
+                maskX = E{1}*maskXSmall*E{1}' + E{2}*maskXSmall*E{2}';
             end
             for r = 1:m(2)
                 p = ind(:,r);
@@ -254,8 +254,14 @@
                 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}';
+                maskYSmall = kron(speye(m(1), m(1)), mask);
+
+                maskY = E{1}*maskYSmall*E{1}' + E{2}*maskYSmall*E{2}';
+                mask = maskX + maskY;
+                mask = mask>0;
+
+                maskSmall = maskXSmall + maskYSmall;
+                maskSmall = maskSmall>0;
             end
             for r = 1:m(1)
                 p = ind(r,:);
@@ -292,7 +298,11 @@
                                 D = D + E{j}*D2_temp{j,k,l}*E{l}';
                                 D2_temp{j,k,l} = [];
                             else
-                                D = D + E{j}*D1{i}*C_mat{i,j,k,l}*D1{k}*E{l}';
+                                if hollow
+                                    D = D + E{j}*(maskSmall*D1{i})*C_mat{i,j,k,l}*D1{k}*E{l}';
+                                else
+                                    D = D + E{j}*(D1{i})*C_mat{i,j,k,l}*D1{k}*E{l}';
+                                end
                             end
                         end
                     end