Mercurial > repos > public > sbplib

diff +scheme/Schrodinger2dCurve.m @ 694:1157375c678a feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add oposit jacobian to penalty term
author Ylva Rydin <ylva.rydin@telia.com>
date Tue, 19 Sep 2017 09:47:14 +0200
parents f235284e2eb1
children ba0d31ce4121
line wrap: on
line diff
--- a/+scheme/Schrodinger2dCurve.m	Fri Sep 15 09:24:21 2017 +0200
+++ b/+scheme/Schrodinger2dCurve.m	Tue Sep 19 09:47:14 2017 +0200
@@ -100,7 +100,7 @@
             obj.e_n  = kr(I_u,e_r_v);
             obj.du_w = kr(d1_l_u,I_v);
             obj.dv_w = (obj.e_w'*obj.Dv)';
-            obj.du_e = kr(d1_r_u,I_v);
+            obj.du_e = kr(d1_r_u,I_v);   
             obj.dv_e = (obj.e_e'*obj.Dv)';
             obj.du_s = (obj.e_s'*obj.Du)';
             obj.dv_s = kr(I_u,d1_l_v);
@@ -227,8 +227,8 @@
         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_n_u, d_t_u,  coeff_t_u, coeff_n_u,s_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t,gamm_u] = obj.get_boundary_ops(boundary);
-            [e_v, d_n_v, d_t_v,  coeff_t_v, coeff_n_v s_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t,gamm_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
+            [e_u, d_n_u, d_t_u,  coeff_t_u, coeff_n_u,s_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t,gamm_u,Ji_u] = obj.get_boundary_ops(boundary);
+            [e_v, d_n_v, d_t_v,  coeff_t_v, coeff_n_v s_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t,gamm_v,Ji_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
 
             a_n_u = @(t) spdiag(coeff_n_u(t));
             a_t_u = @(t) spdiag(coeff_t_u(t));
@@ -250,8 +250,8 @@
             penalty_parameter_3 = @(t) halfnorm_inv_u_n * e_u  * (gamm(t) );
 
 
-            closure =@(t) sqrt(obj.Ji)*obj.c^2 * ( penalty_parameter_1(t)*e_u' + penalty_parameter_2(t)*F_u(t)' + penalty_parameter_3(t)*e_u')*sqrt(obj.Ji);
-            penalty =@(t) sqrt(obj.Ji)*obj.c^2 * ( -penalty_parameter_1(t)*e_v' - penalty_parameter_2(t)*F_v(t)' - penalty_parameter_3(t)*e_v')*sqrt(obj.Ji);
+            closure =@(t) sqrt(Ji_u)*obj.c^2 * ( penalty_parameter_1(t)*e_u' + penalty_parameter_2(t)*F_u(t)' + penalty_parameter_3(t)*e_u')*sqrt(Ji_u);
+            penalty =@(t) sqrt(Ji_v)*obj.c^2 * ( -penalty_parameter_1(t)*e_v' - penalty_parameter_2(t)*F_v(t)' - penalty_parameter_3(t)*e_v')*sqrt(Ji_v);
         end
 
 
@@ -318,6 +318,7 @@
                     halfnorm_inv_t = obj.Hiu;
                     halfnorm_t = obj.Hu;
             end
+             Ji = obj.Ji;
         end
 
         function N = size(obj)