Mercurial > repos > public > sbplib

diff +scheme/Schrodinger2dCurve.m @ 492:6b95a894cbd7 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fixed a bug in the metric coefficients but something is wrong at the boundaries
author Ylva Rydin <ylva.rydin@telia.com>
date Wed, 15 Feb 2017 11:21:04 +0100
parents 26125cfefe11
children 6b8297f66c91
line wrap: on
line diff
--- a/+scheme/Schrodinger2dCurve.m	Fri Feb 10 14:29:53 2017 +0100
+++ b/+scheme/Schrodinger2dCurve.m	Wed Feb 15 11:21:04 2017 +0100
@@ -19,7 +19,7 @@
         D1_v, D1_u
         D2_v, D2_u
         Du, Dv
-        
+        x,y
  
         e_w, e_e, e_s, e_n
         du_w, dv_w
@@ -125,6 +125,8 @@
             [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2));
             x = reshape(obj.xm,obj.m_tot,1);
             y = reshape(obj.ym,obj.m_tot,1);
+            obj.x=x;
+            obj.y=y;
 
             x_tau = reshape(x_tau,obj.m_tot,1);
             y_tau = reshape(y_tau,obj.m_tot,1); 
@@ -167,11 +169,18 @@
          
             Ji = spdiags(1./J, 0, obj.m_tot, obj.m_tot);
             obj.Ji=Ji;
-            obj.g_1 = x_tau.*y_v-y_tau.*x_v;
-            obj.g_2 = -x_tau.*y_u + y_tau.*x_u;
+            obj.g_1 = x_v.*y_tau-x_tau.*y_v;
+            obj.g_2 = x_tau.*y_u - y_tau.*x_u;
+            
+            b1 = spdiags(obj.g_1, 0, obj.m_tot, obj.m_tot);
+            b2 = spdiags(obj.g_2, 0, obj.m_tot, obj.m_tot);
+            
+            b1_u = spdiags(obj.Du*obj.g_1, 0, obj.m_tot, obj.m_tot);
+            b2_v = spdiags(obj.Dv*obj.g_2, 0, obj.m_tot, obj.m_tot);
             
             %Add the flux splitting
-            D = Ji*(obj.g_1.*obj.Du + obj.g_2.*obj.Dv + 1i*obj.c^2*(Duu + Duv + Dvu + Dvv)); %% g_1' och g_2'?
+          % D = Ji*(-b1*obj.Du  -b2*obj.Dv + 1i*obj.c^2*(Duu + Duv + Dvu + Dvv)); 
+             D = Ji*(-1/2*(b1*obj.Du-b1_u+obj.Du*b1) - 1/2*(b2*obj.Dv - b2_v +obj.Dv*b2) + 1i*obj.c^2*(Duu + Duv + Dvu + Dvv)); 
             
 %             obj.gamm_u = h_u*ops_u.borrowing.M.d1;
 %             obj.gamm_v = h_v*ops_v.borrowing.M.d1;
@@ -194,12 +203,11 @@
                     a = spdiag(g);
                     tau2 =  (-1*s*a - abs(a))/4;
 
-                    penalty_parameter_1 = 1*1i*halfnorm_inv_n*halfnorm_inv_t*e*F'*halfnorm_t*e;
-                    penalty_parameter_2 = halfnorm_inv_n*e*tau2;
+                    penalty_parameter_1 = 1*1i*halfnorm_inv_n*halfnorm_inv_t*F*e'*halfnorm_t*e;
+                     penalty_parameter_2 = halfnorm_inv_n*e*tau2;
 
-                    closure = obj.Ji*obj.c^2 * penalty_parameter_1*e' +obj.Ji* penalty_parameter_2*e';
-                   % penalty = -obj.Ji*obj.c^2 * penalty_parameter_2;
-                   penalty=0;
+                    closure = obj.Ji*obj.c^2 * penalty_parameter_1*e' + obj.Ji* penalty_parameter_2*e';
+                    penalty = -obj.Ji*obj.c^2 * penalty_parameter_1*e'+  obj.Ji*penalty_parameter_2*e';
                 
         end