diff +scheme/Schrodinger1dCurve.m @ 495:b91d23271481 feature/quantumTriangles

Added new penalty parameter to the interface in shrodinger curve
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 24 Feb 2017 09:04:02 +0100
parents 6b8297f66c91
children 437fad4a47e1
line wrap: on
line diff
--- a/+scheme/Schrodinger1dCurve.m	Thu Feb 23 09:32:30 2017 +0100
+++ b/+scheme/Schrodinger1dCurve.m	Fri Feb 24 09:04:02 2017 +0100
@@ -24,7 +24,7 @@
     
     methods
         % Solving SE in the form u_t = i*u_xx +i*V on deforming 1D domain;
-        function obj = Schrodinger1dCurve(m,order,V,constJi)
+        function obj = Schrodinger1dCurve(g,m,order,V,constJi)
             default_arg('V',0);
             default_arg('constJi',false)
             xilim={0 1};
@@ -45,6 +45,7 @@
             obj.e_r = ops.e_r;
             obj.d1_l = ops.d1_l;
             obj.d1_r = ops.d1_r;
+            obj.grid = g;
             
             
             if isa(V,'function_handle')
@@ -107,17 +108,18 @@
         function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
             % u denotes the solution in the own domain
             % v denotes the solution in the neighbour domain
-                         [e_u,d_u,s_u] = obj.get_boundary_ops(boundary);
-                         [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
-            
-                         a =  -s_u* 1/2 * 1i ;
-                         b =  a';
+                         [e_u,d_u,s_u,p_u] = obj.get_boundary_ops(boundary);
+                         [e_v,d_v,s_v,p_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
             
-                         tau = b*d_u;
-                         sig = -a*e_u;
+                         a1 =  -s_u* 1/2 * 1i ;
+                         b1 =  a1';
+                         gamma = -1*s*a(p_u,p_u)/2;
             
-                         closure = obj.Hi * (tau*e_u' + sig*d_u');
-                         penalty = obj.Hi * (-tau*e_v' - sig*d_v');
+                         tau = b1*d_u;
+                         sig = -a1*e_u;
+            
+                         closure = @(a) obj.Hi * (tau*e_u' + sig*d_u' + gamma(a));
+                         penalty = @(a) obj.Hi * (-tau*e_v' - sig*d_v' - gamma(a));
         end
         
         % Ruturns the boundary ops and sign for the boundary specified by the string boundary.