--- a/+scheme/ViscoElastic2d.m	Wed Jul 22 15:35:25 2020 -0700
+++ b/+scheme/ViscoElastic2d.m	Sat Jul 25 15:56:31 2020 -0700
@@ -42,6 +42,11 @@
         n_w, n_e, n_s, n_n % Physical normals
         tangent_w, tangent_e, tangent_s, tangent_n % Physical tangents
 
+        tau1_w, tau1_e, tau1_s, tau1_n  % Return scalar field
+        tau2_w, tau2_e, tau2_s, tau2_n  % Return scalar field
+        tau_n_w, tau_n_e, tau_n_s, tau_n_n % Return scalar field
+        tau_t_w, tau_t_e, tau_t_s, tau_t_n % Return scalar field
+
         elasticObj
 
     end
@@ -224,6 +229,57 @@
             obj.e_scalar_s = elasticObj.e_scalar_s;
             obj.e_scalar_n = elasticObj.e_scalar_n;
 
+            % -- Create new traction operators including viscous strain contribution --
+            tau1 = struct;
+            tau2 = struct;
+            tau_n = struct;
+            tau_t = struct;
+            boundaries = {'w', 'e', 's', 'n'};
+            for bNumber = 1:numel(boundaries)
+                b = boundaries{bNumber};
+
+                n = elasticObj.getNormal(b);
+                t = elasticObj.getTangent(b);
+                e = elasticObj.getBoundaryOperatorForScalarField('e', b);
+                tau1.(b) = (elasticObj.getBoundaryOperator('tau1', b)'*Eu')';
+                tau2.(b) = (elasticObj.getBoundaryOperator('tau2', b)'*Eu')';
+
+                % Add viscous contributions
+                for i = 1:dim
+                    for k = 1:dim
+                        for l = 1:dim
+                            tau1.(b) = tau1.(b) - (n{i}*e'*C{i,1,k,l}*eGamma{k,l}')';
+                            tau2.(b) = tau2.(b) - (n{i}*e'*C{i,2,k,l}*eGamma{k,l}')';
+                        end
+                    end
+                end
+
+                % Compute normal and tangential components
+                tau_n.(b) = tau1.(b)*n{1} + tau2.(b)*n{2};
+                tau_t.(b) = tau1.(b)*t{1} + tau2.(b)*t{2};
+            end
+
+            obj.tau1_w = tau1.w;
+            obj.tau1_e = tau1.e;
+            obj.tau1_s = tau1.s;
+            obj.tau1_n = tau1.n;
+
+            obj.tau2_w = tau2.w;
+            obj.tau2_e = tau2.e;
+            obj.tau2_s = tau2.s;
+            obj.tau2_n = tau2.n;
+
+            obj.tau_n_w = tau_n.w;
+            obj.tau_n_e = tau_n.e;
+            obj.tau_n_s = tau_n.s;
+            obj.tau_n_n = tau_n.n;
+
+            obj.tau_t_w = tau_t.w;
+            obj.tau_t_e = tau_t.e;
+            obj.tau_t_s = tau_t.s;
+            obj.tau_t_n = tau_t.n;
+            %----------------------------------------
+
             % Misc.
             obj.elasticObj = elasticObj;
             obj.m = elasticObj.m;
@@ -280,28 +336,26 @@
             switch component
             case 't'
                 dir = obj.getTangent(boundary);
+                tau = obj.getBoundaryOperator('tau_t', boundary);
             case 'n'
                 dir = obj.getNormal(boundary);
+                tau = obj.getBoundaryOperator('tau_n', boundary);
             case 1
                 dir = {1, 0};
+                tau = obj.getBoundaryOperator('tau1', boundary);
             case 2
                 dir = {0, 1};
+                tau = obj.getBoundaryOperator('tau2', boundary);
             end
 
             switch type
             case {'F','f','Free','free','traction','Traction','t','T'}
 
-                % Add viscous part of closure
-                for i = 1:dim
-                    for j = 1:dim
-                        for k = 1:dim
-                            for l = 1:dim
-                                for m = 1:dim
-                                    closure = closure + eU{m}*( (RHO*H)\(C{i,j,k,l}*e*dir{m}*dir{j}*H_gamma*n{i}*e'*eGamma{k,l}') );
-                                end
-                            end
-                        end
-                    end
+                % Set elastic closure to zero
+                closure = 0*closure;
+
+                for m = 1:dim
+                    closure = closure - eU{m}*( (RHO*H)\(e*dir{m}*H_gamma*tau') );
                 end
 
             case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
@@ -514,7 +568,7 @@
             assertIsMember(boundary, {'w', 'e', 's', 'n'})
             assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n', 'tau_t'})
 
-            o = obj.elasticObj.([op, '_', boundary]);
+            o = obj.([op, '_', boundary]);
 
         end