Mercurial > repos > public > sbplib

changeset 1118:07d0caf915e4 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
author Martin Almquist <malmquist@stanford.edu>
date Sun, 05 May 2019 19:05:31 -0700
parents 27019aca2f17
children 8984b12feba6 b8bd54332763
files +scheme/Elastic2dVariable.m
diffstat 1 files changed, 39 insertions(+), 33 deletions(-) [+]
line wrap: on
line diff
--- a/+scheme/Elastic2dVariable.m	Sat May 04 14:44:34 2019 -0700
+++ b/+scheme/Elastic2dVariable.m	Sun May 05 19:05:31 2019 -0700
@@ -61,11 +61,13 @@
     methods
 
         % The coefficients can either be function handles or grid functions
-        function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet)
+        % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
+        function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet, optFlag)
             default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
             default_arg('lambda', @(x,y) 0*x+1);
             default_arg('mu', @(x,y) 0*x+1);
             default_arg('rho', @(x,y) 0*x+1);
+            default_arg('optFlag', false);
             dim = 2;
 
             assert(isa(g, 'grid.Cartesian'))
@@ -356,42 +358,46 @@
             obj.dim = dim;
 
             % B, used for adjoint optimization
-            B = cell(dim, 1);
-            for i = 1:dim
-                B{i} = cell(m_tot, 1);
-            end
+            B = [];
+            if optFlag
+                B = cell(dim, 1);
+                for i = 1:dim
+                    B{i} = cell(m_tot, 1);
+                end
+
+                B0 = sparse(m_tot, m_tot);
+                for i = 1:dim
+                    for j = 1:m_tot
+                        B{i}{j} = B0;
+                    end
+                end
+
+                ind = grid.funcToMatrix(g, 1:m_tot);
 
-            for i = 1:dim
-                for j = 1:m_tot
-                    B{i}{j} = sparse(m_tot, m_tot);
+                % Direction 1
+                for k = 1:m(1)
+                    c = sparse(m(1),1);
+                    c(k) = 1;
+                    [~, B_1D] = ops{1}.D2(c);
+                    for l = 1:m(2)
+                        p = ind(:,l);
+                        B{1}{(k-1)*m(2) + l}(p, p) = B_1D;
+                    end
+                end
+
+                % Direction 2
+                for k = 1:m(2)
+                    c = sparse(m(2),1);
+                    c(k) = 1;
+                    [~, B_1D] = ops{2}.D2(c);
+                    for l = 1:m(1)
+                        p = ind(l,:);
+                        B{2}{(l-1)*m(2) + k}(p, p) = B_1D;
+                    end
                 end
             end
-
-            ind = grid.funcToMatrix(g, 1:m_tot);
+            obj.B = B;
 
-            % Direction 1
-            for k = 1:m(1)
-                c = sparse(m(1),1);
-                c(k) = 1;
-                [~, B_1D] = ops{1}.D2(c);
-                for l = 1:m(2)
-                    p = ind(:,l);
-                    B{1}{(k-1)*m(2) + l}(p, p) = B_1D;
-                end
-            end
-
-            % Direction 2
-            for k = 1:m(2)
-                c = sparse(m(2),1);
-                c(k) = 1;
-                [~, B_1D] = ops{2}.D2(c);
-                for l = 1:m(1)
-                    p = ind(l,:);
-                    B{2}{(l-1)*m(2) + k}(p, p) = B_1D;
-                end
-            end
-
-            obj.B = B;
 
         end