--- a/+scheme/Elastic2dCurvilinearAnisotropic.m	Mon Oct 14 13:54:50 2019 -0700
+++ b/+scheme/Elastic2dCurvilinearAnisotropic.m	Mon Oct 21 18:08:22 2019 -0700
@@ -24,6 +24,7 @@
         D  % Total operator
 
         Dx, Dy % Physical derivatives
+        sigma % Cell matrix of physical stress operators
         n_w, n_e, n_s, n_n % Physical normals
 
         % Boundary operators in cell format, used for BC
@@ -305,6 +306,23 @@
             obj.en_s = (obj.n_s{1}*obj.e1_s' + obj.n_s{2}*obj.e2_s')';
             obj.en_n = (obj.n_n{1}*obj.e1_n' + obj.n_n{2}*obj.e2_n')';
 
+            % Stress operators
+            sigma = cell(dim, dim);
+            D1 = {obj.Dx, obj.Dy};
+            E = obj.E;
+            N = length(obj.RHO);
+            for i = 1:dim
+                for j = 1:dim
+                    sigma{i,j} = sparse(N,2*N);
+                    for k = 1:dim
+                        for l = 1:dim
+                            sigma{i,j} = sigma{i,j} + obj.C{i,j,k,l}*D1{k}*E{l}';
+                        end
+                    end
+                end
+            end
+            obj.sigma = sigma;
+
             % Misc.
             obj.refObj = refObj;
             obj.m = refObj.m;