sbplib: +scheme/Schrodinger2dCurve.m annotate

annotate +scheme/Schrodinger2dCurve.m @ 698:0122cbe2e6d3 feature/optim

Add hamiltonian

author	Ylva Rydin <ylva.rydin@telia.com>
date	Fri, 13 Oct 2017 09:16:33 +0200
parents	ba0d31ce4121
children	8f1eae1450b2

rev	line source
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	1 classdef Schrodinger2dCurve < scheme.Scheme
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	2 properties
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	3 m % Number of points in each direction, possibly a vector
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	4 h % Grid spacing
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	5
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	6 grid
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	7 xm, ym
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	8
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	9 order % Order accuracy for the approximation
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	10
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	11 D % non-stabalized scheme operator
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	12 M % Derivative norm
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	13 H % Discrete norm
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	14 Hi
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	15 H_u, H_v % Norms in the x and y directions
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	16 Hu,Hv % Kroneckerd norms. 1'Hxv corresponds to integration in the x dir.
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	17 Hi_u, Hi_v
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	18 Hiu, Hiv
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	19 D1_v, D1_u
26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	20 D2_v, D2_u
26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	21 Du, Dv
492 6b95a894cbd7 fixed a bug in the metric coefficients but something is wrong at the boundaries Ylva Rydin <ylva.rydin@telia.com> parents: 491 diff changeset	22 x,y
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	23 b1, b2
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	24 b1_u,b2_v
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	25 DU, DV, DUU, DUV, DVU, DVV
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	26
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	27 e_w, e_e, e_s, e_n
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	28 du_w, dv_w
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	29 du_e, dv_e
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	30 du_s, dv_s
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	31 du_n, dv_n
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	32 g_1, g_2
26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	33 c
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	34 ind
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	35 t_up
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	36
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	37 a11, a12, a22
26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	38 m_tot, m_u, m_v
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	39 p
518 4709f2329372 Added J to propertied Ylva Rydin <ylva.rydin@telia.com> parents: 517 diff changeset	40 Ji, J
698 0122cbe2e6d3 Add hamiltonian Ylva Rydin <ylva.rydin@telia.com> parents: 695 diff changeset	41 Hamiltonian
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	42 end
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	43
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	44 methods
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	45 function obj = Schrodinger2dCurve(g ,order, opSet,p)
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	46 default_arg('opSet',@sbp.D2Variable);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	47 default_arg('c', 1);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	48
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	49 obj.p=p;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	50 obj.c=1;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	51
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	52 m = g.size();
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	53 obj.m_u = m(1);
26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	54 obj.m_v = m(2);
26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	55 obj.m_tot = g.N();
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	56 obj.grid=g;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	57
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	58 h = g.scaling();
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	59 h_u = h(1);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	60 h_v = h(2);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	61
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	62 % Operators
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	63 ops_u = opSet(obj.m_u, {0, 1}, order);
26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	64 ops_v = opSet(obj.m_v, {0, 1}, order);
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	65
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	66 I_u = speye(obj.m_u);
26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	67 I_v = speye(obj.m_v);
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	68
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	69 obj.D1_u = ops_u.D1;
26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	70 obj.D2_u = ops_u.D2;
698 0122cbe2e6d3 Add hamiltonian Ylva Rydin <ylva.rydin@telia.com> parents: 695 diff changeset	71
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	72 H_u = ops_u.H;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	73 Hi_u = ops_u.HI;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	74 e_l_u = ops_u.e_l;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	75 e_r_u = ops_u.e_r;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	76 d1_l_u = ops_u.d1_l;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	77 d1_r_u = ops_u.d1_r;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	78
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	79 obj.D1_v = ops_v.D1;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	80 obj.D2_v = ops_v.D2;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	81 H_v = ops_v.H;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	82 Hi_v = ops_v.HI;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	83 e_l_v = ops_v.e_l;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	84 e_r_v = ops_v.e_r;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	85 d1_l_v = ops_v.d1_l;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	86 d1_r_v = ops_v.d1_r;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	87
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	88 obj.Du = kr(obj.D1_u,I_v);
26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	89 obj.Dv = kr(I_u,obj.D1_v);
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	90
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	91 obj.H = kr(H_u,H_v);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	92 obj.Hi = kr(Hi_u,Hi_v);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	93 obj.Hu = kr(H_u,I_v);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	94 obj.Hv = kr(I_u,H_v);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	95 obj.Hiu = kr(Hi_u,I_v);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	96 obj.Hiv = kr(I_u,Hi_v);
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	97
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	98 obj.e_w = kr(e_l_u,I_v);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	99 obj.e_e = kr(e_r_u,I_v);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	100 obj.e_s = kr(I_u,e_l_v);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	101 obj.e_n = kr(I_u,e_r_v);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	102 obj.du_w = kr(d1_l_u,I_v);
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	103 obj.dv_w = (obj.e_w'*obj.Dv)';
694 1157375c678a Add oposit jacobian to penalty term Ylva Rydin <ylva.rydin@telia.com> parents: 520 diff changeset	104 obj.du_e = kr(d1_r_u,I_v);
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	105 obj.dv_e = (obj.e_e'*obj.Dv)';
26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	106 obj.du_s = (obj.e_s'*obj.Du)';
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	107 obj.dv_s = kr(I_u,d1_l_v);
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	108 obj.du_n = (obj.e_n'*obj.Du)';
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	109 obj.dv_n = kr(I_u,d1_r_v);
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	110
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	111 obj.DUU = sparse(obj.m_tot);
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	112 obj.DVV = sparse(obj.m_tot);
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	113 obj.ind = grid.funcToMatrix(obj.grid, 1:obj.m_tot);
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	114
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	115 obj.m = m;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	116 obj.h = [h_u h_v];
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	117 obj.order = order;
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	118 obj.D = @(t)obj.d_fun(t);
698 0122cbe2e6d3 Add hamiltonian Ylva Rydin <ylva.rydin@telia.com> parents: 695 diff changeset	119 obj.Hamiltonian = @(t)obj.h_fun(t);
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	120 obj.variable_update(0);
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	121 end
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	122
698 0122cbe2e6d3 Add hamiltonian Ylva Rydin <ylva.rydin@telia.com> parents: 695 diff changeset	123 function [D,n] = d_fun(obj,t)
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	124 % obj.update_vairables(t); In driscretization?
515 ba92b27da5a0 Tried to make skewsym. Ylva Rydin <ylva.rydin@telia.com> parents: 500 diff changeset	125 % D = obj.Ji(-1/2(obj.b1obj.Du-obj.b1_u+obj.Duobj.b1) -
ba92b27da5a0 Tried to make skewsym. Ylva Rydin <ylva.rydin@telia.com> parents: 500 diff changeset	126 % 1/2(obj.b2obj.Dv - obj.b2_v +obj.Dv*obj.b2) +
ba92b27da5a0 Tried to make skewsym. Ylva Rydin <ylva.rydin@telia.com> parents: 500 diff changeset	127 % 1iobj.c^2(obj.DUU + obj.DUV + obj.DVU + obj.DVV)); (ols
ba92b27da5a0 Tried to make skewsym. Ylva Rydin <ylva.rydin@telia.com> parents: 500 diff changeset	128 % not skew sym disc
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	129
520 f235284e2eb1 changed sign in penalty parameter Ylva Rydin <ylva.rydin@telia.com> parents: 519 diff changeset	130 D = sqrt(obj.Ji)(-1/2(obj.b1obj.Du + obj.Duobj.b1) - 1/2(obj.b2obj.Dv + obj.Dvobj.b2) + 1iobj.c^2(obj.DUU + obj.DUV + obj.DVU + obj.DVV))sqrt(obj.Ji);
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	131 end
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	132
698 0122cbe2e6d3 Add hamiltonian Ylva Rydin <ylva.rydin@telia.com> parents: 695 diff changeset	133 function [Hamiltonian,n] = h_fun(obj,t)
0122cbe2e6d3 Add hamiltonian Ylva Rydin <ylva.rydin@telia.com> parents: 695 diff changeset	134 Hamiltonian = obj.c^2(obj.DUU + obj.DUV + obj.DVU + obj.DVV)sqrt(obj.Ji);
0122cbe2e6d3 Add hamiltonian Ylva Rydin <ylva.rydin@telia.com> parents: 695 diff changeset	135 end
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	136
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	137 function [D ]= variable_update(obj,t)
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	138 % Metric derivatives
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	139 if(obj.t_up == t)
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	140 return
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	141 else
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	142 ti = parametrization.Ti.points(obj.p{1}(t),obj.p{2}(t),obj.p{3}(t),obj.p{4}(t));
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	143 ti_tau = parametrization.Ti.points(obj.p{5}(t),obj.p{6}(t),obj.p{7}(t),obj.p{8}(t));
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	144 lcoords=points(obj.grid);
498 324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	145 [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_v,2));
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	146 [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_v,2));
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	147 x = reshape(obj.xm,obj.m_tot,1);
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	148 y = reshape(obj.ym,obj.m_tot,1);
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	149 obj.x = x;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	150 obj.y = y;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	151
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	152 x_tau = reshape(x_tau,obj.m_tot,1);
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	153 y_tau = reshape(y_tau,obj.m_tot,1);
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	154
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	155 x_u = obj.Du*x;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	156 x_v = obj.Dv*x;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	157 y_u = obj.Du*y;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	158 y_v = obj.Dv*y;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	159
500 83734c26b8e3 Small changes Ylva Rydin <ylva.rydin@telia.com> parents: 499 diff changeset	160 J = (x_u.y_v - x_v.y_u);
518 4709f2329372 Added J to propertied Ylva Rydin <ylva.rydin@telia.com> parents: 517 diff changeset	161 obj.J = spdiags(J, 0, obj.m_tot, obj.m_tot);
500 83734c26b8e3 Small changes Ylva Rydin <ylva.rydin@telia.com> parents: 499 diff changeset	162
83734c26b8e3 Small changes Ylva Rydin <ylva.rydin@telia.com> parents: 499 diff changeset	163 Ji = spdiags(1./J, 0, obj.m_tot, obj.m_tot);
83734c26b8e3 Small changes Ylva Rydin <ylva.rydin@telia.com> parents: 499 diff changeset	164 obj.Ji = Ji;
83734c26b8e3 Small changes Ylva Rydin <ylva.rydin@telia.com> parents: 499 diff changeset	165
83734c26b8e3 Small changes Ylva Rydin <ylva.rydin@telia.com> parents: 499 diff changeset	166 a11 = Ji* (x_v.^2 + y_v.^2);
83734c26b8e3 Small changes Ylva Rydin <ylva.rydin@telia.com> parents: 499 diff changeset	167 a12 = -Ji* (x_u.x_v + y_u.y_v);
83734c26b8e3 Small changes Ylva Rydin <ylva.rydin@telia.com> parents: 499 diff changeset	168 a22 = Ji* (x_u.^2 + y_u.^2);
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	169
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	170 obj.a11 = a11;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	171 obj.a12 = a12;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	172 obj.a22 = a22;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	173
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	174 % Assemble full operators
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	175 L_12 = spdiags(a12, 0, obj.m_tot, obj.m_tot);
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	176 obj.DUV = obj.DuL_12obj.Dv;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	177 obj.DVU = obj.DvL_12obj.Du;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	178
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	179
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	180 for i = 1:obj.m_v
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	181 D = obj.D2_u(a11(obj.ind(:,i)));
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	182 p = obj.ind(:,i);
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	183 obj.DUU(p,p) = D;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	184 end
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	185
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	186 for i = 1:obj.m_u
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	187 D = obj.D2_v(a22(obj.ind(i,:)));
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	188 p = obj.ind(i,:);
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	189 obj.DVV(p,p) = D;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	190 end
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	191
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	192 obj.g_1 = x_v.y_tau-x_tau.y_v;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	193 obj.g_2 = x_tau.y_u - y_tau.x_u;
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	194
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	195 obj.b1 = spdiags(obj.g_1, 0, obj.m_tot, obj.m_tot);
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	196 obj.b2 = spdiags(obj.g_2, 0, obj.m_tot, obj.m_tot);
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	197
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	198 obj.b1_u = spdiags(obj.Du*obj.g_1, 0, obj.m_tot, obj.m_tot);
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	199 obj.b2_v = spdiags(obj.Dv*obj.g_2, 0, obj.m_tot, obj.m_tot);
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	200 obj.t_up=t;
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	201 end
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	202 end
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	203
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	204 % Closure functions return the opertors applied to the own doamin to close the boundary
516 afff85574ddb Changed the penalty on the first derivative penalty terms, did not solve the problem Ylva Rydin <ylva.rydin@telia.com> parents: 515 diff changeset	205
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	206 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	207 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	208 % type is a string specifying the type of boundary condition if there are several.
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	209 % data is a function returning the data that should be applied at the boundary.
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	210 % neighbour_scheme is an instance of Scheme that should be interfaced to.
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	211 % neighbour_boundary is a string specifying which boundary to interface to.
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	212 function [closure, penalty] = boundary_condition(obj, boundary,~)
499 f1465e6aeb26 Commit before changing branch Ylva Rydin <ylva.rydin@telia.com> parents: 498 diff changeset	213 [e, d_n, d_t, coeff_t, coeff_n s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary);
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	214
498 324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	215 a_t = @(t) spdiag(coeff_t(t));
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	216 a_n = @(t) spdiag(coeff_n(t));
499 f1465e6aeb26 Commit before changing branch Ylva Rydin <ylva.rydin@telia.com> parents: 498 diff changeset	217
498 324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	218
499 f1465e6aeb26 Commit before changing branch Ylva Rydin <ylva.rydin@telia.com> parents: 498 diff changeset	219 F = @(t)(s * a_n(t)d_n' + s a_t(t) *d_t')';
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	220 tau1 = 1;
498 324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	221 a = @(t)spdiag(g(t));
517 7a091a3527df sign change in SAT-TERM Ylva Rydin <ylva.rydin@telia.com> parents: 516 diff changeset	222 tau2 = @(t) (1sa(t))/2;
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	223
498 324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	224 penalty_parameter_1 = @(t) 11ihalfnorm_inv_nhalfnorm_inv_tF(t)e'halfnorm_t*e;
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	225 penalty_parameter_2 = @(t) halfnorm_inv_netau2(t);
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	226
515 ba92b27da5a0 Tried to make skewsym. Ylva Rydin <ylva.rydin@telia.com> parents: 500 diff changeset	227 closure = @(t) sqrt(obj.Ji)(obj.c^2 penalty_parameter_1(t)e' + penalty_parameter_2(t)e')*sqrt(obj.Ji);
520 f235284e2eb1 changed sign in penalty parameter Ylva Rydin <ylva.rydin@telia.com> parents: 519 diff changeset	228 penalty = @(t) -sqrt(obj.Ji)(obj.c^2 penalty_parameter_1(t)e' + penalty_parameter_2(t)e')*sqrt(obj.Ji);
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	229
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	230 end
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	231
26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	232 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	233 % u denotes the solution in the own domain
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	234 % v denotes the solution in the neighbour domain
694 1157375c678a Add oposit jacobian to penalty term Ylva Rydin <ylva.rydin@telia.com> parents: 520 diff changeset	235 [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);
1157375c678a Add oposit jacobian to penalty term Ylva Rydin <ylva.rydin@telia.com> parents: 520 diff changeset	236 [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);
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	237
498 324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	238 a_n_u = @(t) spdiag(coeff_n_u(t));
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	239 a_t_u = @(t) spdiag(coeff_t_u(t));
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	240 a_n_v = @(t) spdiag(coeff_n_v(t));
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	241 a_t_v = @(t) spdiag(coeff_t_v(t));
499 f1465e6aeb26 Commit before changing branch Ylva Rydin <ylva.rydin@telia.com> parents: 498 diff changeset	242
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	243
499 f1465e6aeb26 Commit before changing branch Ylva Rydin <ylva.rydin@telia.com> parents: 498 diff changeset	244 F_u = @(t)(a_n_u(t)d_n_u' + a_t_u(t)d_t_u')';
f1465e6aeb26 Commit before changing branch Ylva Rydin <ylva.rydin@telia.com> parents: 498 diff changeset	245 F_v = @(t)(a_n_v(t)d_n_v' + a_t_v(t)d_t_v')';
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	246
498 324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	247 a = @(t)spdiag(gamm_u(t));
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	248
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	249 tau = s_u11i/2;
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	250 sig = -s_u11i/2;
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	251 gamm = @(t) (-s_u*a(t))/2;
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	252
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	253 penalty_parameter_1 = @(t) halfnorm_inv_u_n(tauhalfnorm_inv_u_tF_u(t)e_u'halfnorm_u_te_u);
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	254 penalty_parameter_2 = @(t) halfnorm_inv_u_n * e_u * (sig );
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	255 penalty_parameter_3 = @(t) halfnorm_inv_u_n * e_u * (gamm(t) );
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	256
4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	257
694 1157375c678a Add oposit jacobian to penalty term Ylva Rydin <ylva.rydin@telia.com> parents: 520 diff changeset	258 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);
695 ba0d31ce4121 Changed jacobian at right side of the penalty term in 2D, seems like it is working Ylva Rydin <ylva.rydin@telia.com> parents: 694 diff changeset	259 penalty =@(t) sqrt(Ji_u)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);
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	260 end
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	261
497 4905446f165e Added 2D interface to shrodinger Ylva Rydin <ylva.rydin@telia.com> parents: 493 diff changeset	262
695 ba0d31ce4121 Changed jacobian at right side of the penalty term in 2D, seems like it is working Ylva Rydin <ylva.rydin@telia.com> parents: 694 diff changeset	263 function [e, d_n, d_t, coeff_t,coeff_n, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g,Ji] = get_boundary_ops(obj, boundary)
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	264
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	265 % gridMatrix = zeros(obj.m(2),obj.m(1));
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	266 % gridMatrix(:) = 1:numel(gridMatrix);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	267
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	268 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	269
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	270 switch boundary
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	271 case 'w'
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	272 e = obj.e_w;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	273 d_n = obj.du_w;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	274 d_t = obj.dv_w;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	275 s = -1;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	276
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	277 I = ind(1,:);
498 324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	278 coeff_t = @(t)obj.a12(I);
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	279 coeff_n =@(t) obj.a11(I);
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	280 g = @(t)obj.g_1(I);
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	281 case 'e'
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	282 e = obj.e_e;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	283 d_n = obj.du_e;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	284 d_t = obj.dv_e;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	285 s = 1;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	286
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	287 I = ind(end,:);
498 324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	288 coeff_t = @(t)obj.a12(I);
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	289 coeff_n = @(t)obj.a11(I);
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	290 g = @(t)obj.g_1(I);
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	291 case 's'
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	292 e = obj.e_s;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	293 d_n = obj.dv_s;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	294 d_t = obj.du_s;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	295 s = -1;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	296
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	297 I = ind(:,1)';
498 324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	298 coeff_t = @(t)obj.a12(I);
499 f1465e6aeb26 Commit before changing branch Ylva Rydin <ylva.rydin@telia.com> parents: 498 diff changeset	299 coeff_n = @(t)obj.a22(I);
498 324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	300 g = @(t)obj.g_2(I);
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	301 case 'n'
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	302 e = obj.e_n;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	303 d_n = obj.dv_n;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	304 d_t = obj.du_n;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	305 s = 1;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	306
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	307 I = ind(:,end)';
498 324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	308 coeff_t = @(t)obj.a12(I);
499 f1465e6aeb26 Commit before changing branch Ylva Rydin <ylva.rydin@telia.com> parents: 498 diff changeset	309 coeff_n = @(t)obj.a22(I);
498 324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	310 g = @(t)obj.g_2(I);
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	311 otherwise
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	312 error('No such boundary: boundary = %s',boundary);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	313 end
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	314
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	315 switch boundary
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	316 case {'w','e'}
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	317 halfnorm_inv_n = obj.Hiu;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	318 halfnorm_inv_t = obj.Hiv;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	319 halfnorm_t = obj.Hv;
491 26125cfefe11 Schrodinger 2d runs without exploding with fixed coordinates Ylva Rydin <ylva.rydin@telia.com> parents: 490 diff changeset	320
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	321 case {'s','n'}
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	322 halfnorm_inv_n = obj.Hiv;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	323 halfnorm_inv_t = obj.Hiu;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	324 halfnorm_t = obj.Hu;
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	325 end
694 1157375c678a Add oposit jacobian to penalty term Ylva Rydin <ylva.rydin@telia.com> parents: 520 diff changeset	326 Ji = obj.Ji;
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	327 end
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	328
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	329 function N = size(obj)
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	330 N = prod(obj.m);
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	331 end
b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	332 end
498 324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	333 methods(Static)
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	334 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	335 % and bound_v of scheme schm_v.
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	336 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	337 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	338 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	339 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	340 end
324c927d8b1d chnaged sbp interfacein 1d among many things Ylva Rydin <ylva.rydin@telia.com> parents: 497 diff changeset	341 end
490 b13d44271ead Schrodinger2dCurve Added Ylva Rydin <ylva.rydin@telia.com> parents: diff changeset	342 end

Mercurial > repos > public > sbplib

annotate +scheme/Schrodinger2dCurve.m @ 698:0122cbe2e6d3 feature/optim