Mercurial > repos > public > sbplib

changeset 5:5f6b0b6a012b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First try at interface implementation in WaveCurve2s
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 21 Sep 2015 09:42:15 +0200
parents 8e14b5a577a6
children c3c95b7a1a1c
files +scheme/Wave2dCurve.m
diffstat 1 files changed, 15 insertions(+), 17 deletions(-) [+]
line wrap: on
line diff
--- a/+scheme/Wave2dCurve.m	Fri Sep 18 15:12:44 2015 +0200
+++ b/+scheme/Wave2dCurve.m	Mon Sep 21 09:42:15 2015 +0200
@@ -229,36 +229,34 @@
         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
+            tuning = 1.2;
             [e_u, d_n_u, d_t_u, coeff_n_u, coeff_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t] = obj.get_boundary_ops(boundary);
             [e_v, d_n_v, d_t_v, coeff_n_v, coeff_t_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t] = neighbour_scheme.get_boundary_ops(boundary);
 
             F_u = s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u';
             F_v = s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v';
 
-            tau  = 111111111111111111111111111111;
+            u = obj;
+            v = neighbour_scheme;
+
+            b1_u = gamm_u*u.lambda./u.a11.^2;
+            b2_u = gamm_u*u.lambda./u.a22.^2;
+            b1_v = gamm_v*v.lambda./v.a11.^2;
+            b2_v = gamm_v*v.lambda./v.a22.^2;
+
+
+            tau  = -1./(4*b1_u) -1/(4*b1_v) -1/(4*b2_u) -1/(4*b2_v);
+            m_tot = obj.m(1)*obj.m(2);
+            tau = tuning * spdiags(tau(:),0,m_tot,m_tot);
             sig1 = 1/2;
             sig2 = -1/2;
 
             penalty_parameter_1 = s_u*halfnorm_inv_u_n*(tau + sig1*halfnorm_inv_u_t*F_u'*halfnorm_u_t)*e_u;
             penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u;
 
-            tuning = 1.2;
 
-            alpha_u = obj.alpha;
-            alpha_v = neighbour_scheme.alpha;
-
-            % tau1 < -(alpha_u/gamm_u + alpha_v/gamm_v)
-
-            tau1 = -(alpha_u/gamm_u + alpha_v/gamm_v) * tuning;
-            tau2 = s_u*1/2*alpha_u;
-            sig1 = s_u*(-1/2);
-            sig2 = 0;
-
-            tau = tau1*e_u + tau2*d_u;
-            sig = sig1*e_u + sig2*d_u;
-
-            closure = halfnorm_inv*( tau*e_u' + sig*alpha_u*d_u');
-            penalty = halfnorm_inv*(-tau*e_v' - sig*alpha_v*d_v');
+            closure = obj.Ji*obj.c^2 * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u');
+            penalty = obj.Ji*obj.c^2 * (-penalty_parameter_1*e_v' - penalty_parameter_2*F_v');
         end
 
         % Ruturns the boundary ops and sign for the boundary specified by the string boundary.