Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2dCurve.m @ 369:9d1fc984f40d feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added some comments and cleaned up the code a little
author Ylva Rydin <ylva.rydin@telia.com>
date Thu, 26 Jan 2017 09:57:24 +0100
parents 9b3d7fc61a36
children 459eeb99130f
comparison
equal deleted inserted replaced
368:53abf04f5e4e 369:9d1fc984f40d
1 classdef Hypsyst2dCurve < scheme.Scheme 1 classdef Hypsyst2dCurve < scheme.Scheme
2 properties 2 properties
3 m % Number of points in each direction, possibly a vector 3 m % Number of points in each direction, possibly a vector
4 n %size of system 4 n % size of system
5 h % Grid spacing 5 h % Grid spacing
6 X,Y % Values of x and y for each grid point 6 X,Y % Values of x and y for each grid point
7 7
8 J, Ji %Jacobaian and inverse Jacobian 8 J, Ji % Jacobaian and inverse Jacobian
9 xi,eta 9 xi,eta
10 Xi,Eta 10 Xi,Eta
11 11
12 A,B 12 A,B
13 X_eta, Y_eta 13 X_eta, Y_eta
14 X_xi,Y_xi 14 X_xi,Y_xi
15 order % Order accuracy for the approximation 15 order % Order accuracy for the approximation
16 16
17 D % non-stabalized scheme operator 17 D % non-stabalized scheme operator
18 Ahat, Bhat, E 18 Ahat, Bhat, E
19 19
20 H % Discrete norm 20 H % Discrete norm
21 Hxii,Hetai % Kroneckerd norms in xi and eta. 21 Hxii,Hetai % Kroneckerd norms in xi and eta.
22 I_xi,I_eta, I_N, onesN 22 I_xi,I_eta, I_N, onesN
23 e_w, e_e, e_s, e_n 23 e_w, e_e, e_s, e_n
24 index_w, index_e,index_s,index_n 24 index_w, index_e,index_s,index_n
25 params %parameters for the coeficient matrice 25 params % Parameters for the coeficient matrice
26 end 26 end
27 27
28 28
29 methods 29 methods
30 % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu
30 function obj = Hypsyst2dCurve(m, order, A, B, E, params, ti) 31 function obj = Hypsyst2dCurve(m, order, A, B, E, params, ti)
31 default_arg('E', []) 32 default_arg('E', [])
32 xilim = {0 1}; 33 xilim = {0 1};
33 etalim = {0 1}; 34 etalim = {0 1};
34 35
35 if length(m) == 1 36 if length(m) == 1
36 m = [m m]; 37 m = [m m];
37 end 38 end
38 obj.params = params; 39 obj.params = params;
39 obj.A=A; 40 obj.A=A;
40 obj.B=B; 41 obj.B=B;
41 42
42 obj.Ahat=@(params,x,y,x_eta,y_eta)(A(params,x,y).*y_eta-B(params,x,y).*x_eta); 43 obj.Ahat=@(params,x,y,x_eta,y_eta)(A(params,x,y).*y_eta-B(params,x,y).*x_eta);
43 obj.Bhat=@(params,x,y,x_xi,y_xi)(B(params,x,y).*x_xi-A(params,x,y).*y_xi); 44 obj.Bhat=@(params,x,y,x_xi,y_xi)(B(params,x,y).*x_xi-A(params,x,y).*y_xi);
44 obj.E=@(params,x,y,~,~)E(params,x,y); 45 obj.E=@(params,x,y,~,~)E(params,x,y);
45 46
46 m_xi = m(1); 47 m_xi = m(1);
47 m_eta = m(2); 48 m_eta = m(2);
48 m_tot=m_xi*m_eta; 49 m_tot=m_xi*m_eta;
49 50
50 ops_xi = sbp.D2Standard(m_xi,xilim,order); 51 ops_xi = sbp.D2Standard(m_xi,xilim,order);
51 ops_eta = sbp.D2Standard(m_eta,etalim,order); 52 ops_eta = sbp.D2Standard(m_eta,etalim,order);
52 53
53 obj.xi = ops_xi.x; 54 obj.xi = ops_xi.x;
54 obj.eta = ops_eta.x; 55 obj.eta = ops_eta.x;
55 56
56 obj.Xi = kr(obj.xi,ones(m_eta,1)); 57 obj.Xi = kr(obj.xi,ones(m_eta,1));
57 obj.Eta = kr(ones(m_xi,1),obj.eta); 58 obj.Eta = kr(ones(m_xi,1),obj.eta);
58 59
59 obj.n = length(A(obj.params,0,0)); 60 obj.n = length(A(obj.params,0,0));
60 obj.onesN=ones(obj.n); 61 obj.onesN=ones(obj.n);
61 62
62 obj.index_w=1:m_eta; 63 obj.index_w=1:m_eta;
63 obj.index_e=(m_tot-m_eta+1):m_tot; 64 obj.index_e=(m_tot-m_e
65
66 metric_termsta+1):m_tot;
64 obj.index_s=1:m_eta:(m_tot-m_eta+1); 67 obj.index_s=1:m_eta:(m_tot-m_eta+1);
65 obj.index_n=(m_eta):m_eta:m_tot; 68 obj.index_n=(m_eta):m_eta:m_tot;
66 69
67 I_n = eye(obj.n); 70 I_n = eye(obj.n);
68 I_xi = speye(m_xi); 71 I_xi = speye(m_xi);
69 obj.I_xi = I_xi; 72 obj.I_xi = I_xi;
70 I_eta = speye(m_eta); 73 I_eta = speye(m_eta);
71 obj.I_eta = I_eta; 74 obj.I_eta = I_eta;
72 75
73 D1_xi = kr(I_n, ops_xi.D1, I_eta); 76 D1_xi = kr(I_n, ops_xi.D1, I_eta);
74 obj.Hxii = kr(I_n, ops_xi.HI, I_eta); 77 obj.Hxii = kr(I_n, ops_xi.HI, I_eta);
75 D1_eta = kr(I_n, I_xi, ops_eta.D1); 78 D1_eta = kr(I_n, I_xi, ops_eta.D1);
76 obj.Hetai = kr(I_n, I_xi, ops_eta.HI); 79 obj.Hetai = kr(I_n, I_xi, ops_eta.HI);
77 80
78 obj.e_w = kr(I_n, ops_xi.e_l, I_eta); 81 obj.e_w = kr(I_n, ops_xi.e_l, I_eta);
79 obj.e_e = kr(I_n, ops_xi.e_r, I_eta); 82 obj.e_e = kr(I_n, ops_xi.e_r, I_eta);
80 obj.e_s = kr(I_n, I_xi, ops_eta.e_l); 83 obj.e_s = kr(I_n, I_xi, ops_eta.e_l);
81 obj.e_n = kr(I_n, I_xi, ops_eta.e_r); 84 obj.e_n = kr(I_n, I_xi,
85
86 metric_termsops_eta.e_r);
82 87
83 [X,Y] = ti.map(obj.xi,obj.eta); 88 [X,Y] = ti.map(obj.xi,obj.eta);
84 89
85 [x_xi,x_eta] = gridDerivatives(X,ops_xi.D1,ops_eta.D1); 90 [x_xi,x_eta] = gridDerivatives(X,ops_xi.D1,ops_eta.D1);
86 [y_xi,y_eta] = gridDerivatives(Y,ops_xi.D1, ops_eta.D1); 91 [y_xi,y_eta] = gridDerivatives(Y,ops_xi.D1, ops_eta.D1);
87 92
88 obj.X=reshape(X,m_tot,1); 93 obj.X = reshape(X,m_tot,1);
89 obj.Y=reshape(Y,m_tot,1); 94 obj.Y = reshape(Y,m_tot,1);
90 obj.X_xi=reshape(x_xi,m_tot,1); 95 obj.X_xi = reshape(x_xi,m_tot,1);
91 obj.Y_xi=reshape(y_xi,m_tot,1); 96 obj.Y_xi = reshape(y_xi,m_tot,1);
92 obj.X_eta=reshape(x_eta,m_tot,1); 97 obj.X_eta = reshape(x_eta,m_tot,1);
93 obj.Y_eta=reshape(y_eta,m_tot,1); 98 obj.Y_eta = reshape(y_eta,m_tot,1);
94 99
95 Ahat_evaluated = obj.evaluateCoefficientMatrix(obj.Ahat, obj.X, obj.Y,obj.X_eta,obj.Y_eta); 100 Ahat_evaluated = obj.evaluateCoefficientMatrix(obj.Ahat, obj.X, obj.Y,obj.X_eta,obj.Y_eta);
96 Bhat_evaluated = obj.evaluateCoefficientMatrix(obj.Bhat, obj.X, obj.Y,obj.X_xi,obj.Y_xi); 101 Bhat_evaluated = obj.evaluateCoefficientMatrix(obj.Bhat, obj.X, obj.Y,obj.X_xi,obj.Y_xi);
97 E_evaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,[],[]); 102 E_evaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,[],[]);
98 103
99 obj.m=m; 104 obj.m = m;
100 obj.h=[ops_xi.h ops_eta.h]; 105 obj.h = [ops_xi.h ops_eta.h];
101 obj.order=order; 106 obj.order = order;
102 obj.J=obj.X_xi.*obj.Y_eta - obj.X_eta.*obj.Y_xi; 107 obj.J = obj.X_xi.*obj.Y_eta - obj.X_eta.*obj.Y_xi;
103 obj.Ji =kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot)); 108 obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot));
104 109
105 obj.D=obj.Ji*(-Ahat_evaluated*D1_xi-Bhat_evaluated*D1_eta)-E_evaluated; 110 obj.D = obj.Ji*(-Ahat_evaluated*D1_xi-Bhat_evaluated*D1_eta)-E_evaluated;
106 111
107 end 112 end
108 113
109 % Closure functions return the opertors applied to the own doamin to close the boundary 114 % Closure functions return the opertors applied to the own doamin to close the boundary
110 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 115 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
111 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 116 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w',General boundary conditions'n','s'.
112 % type is a string specifying the type of boundary condition if there are several. 117 % type is a string specifying the type of boundary condition if there are several.
113 % data is a function returning the data that should be applied at the boundary. 118 % data is a function returning the data that should be applied at the boundary.
114 function [closure, penalty] = boundary_condition(obj,boundary,type,L) 119 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
115 default_arg('type','char'); 120 default_arg('type','char');
116 switch type 121 switch type
117 case{'c','char'} 122 case{'c','char'}
118 [closure,penalty]=boundary_condition_char(obj,boundary); 123 [closure,penalty] = boundary_condition_char(obj,boundary);
119 case{'general'} 124 case{'general'}
120 [closure,penalty]=boundary_condition_general(obj,boundary,L); 125 [closure,penalty] = boundary_condition_general(obj,boundary,L);
121 otherwise 126 otherwise
122 error('No such boundary condition') 127 error('No such boundary condition')
123 end 128 end
124 end 129 end
125 130
126 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 131 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundaryGeneral boundary conditions)
127 error('An interface function does not exist yet'); 132 error('An interface function does not exist yet');
128 end 133 end
129 134
130 function N = size(obj) 135 function N = size(obj)
131 N = obj.m; 136 N = obj.m;
132 end 137 end
133 138
134 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_) 139 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_)
135 params=obj.params; 140 params = obj.params;
136 141
137 if isa(mat,'function_handle') 142 if isa(mat,'function_handle')
138 [rows,cols]=size(mat(params,0,0,0,0)); 143 [rows,cols] = size(mat(params,0,0,0,0));
139 x_=kr(obj.onesN,x_); 144 x_ = kr(obj.onesN,x_);
140 y_=kr(obj.onesN,y_); 145 y_ = kr(obj.onesN,y_);
141 matVec=mat(params,X',Y',x_',y_'); 146 matVec = mat(params,X',Y',x_',y_');
142 matVec=sparse(matVec); 147 matVec = sparse(matVec);
143 side=max(length(X),length(Y)); 148 side = max(length(X),length(Y));
144 else 149 else
145 matVec=mat; 150 matVec = mat;
146 [rows,cols]=size(matVec); 151 [rows,cols] = size(matVec);
147 side=max(length(X),length(Y)); 152 side = max(length(X),length(Y));
148 cols=cols/side; 153 cols = cols/side;
149 end 154 end
150 ret=cell(rows,cols); 155
151 156 ret = cell(rows,cols);
152 for ii=1:rows 157 for ii = 1:rows
153 for jj=1:cols 158 for jj = 1:cols
154 ret{ii,jj}=diag(matVec(ii,(jj-1)*side+1:jj*side)); 159 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side));
155 end 160 end
156 end 161 end
157 ret=cell2mat(ret); 162 ret = cell2mat(ret);
158 end 163 end
159 164
160 165 %Characteristic boundary conditions
161 function [closure, penalty]=boundary_condition_char(obj,boundary) 166 function [closure, penalty] = boundary_condition_char(obj,boundary)
162 params=obj.params; 167 params = obj.params;
163 X=obj.X; Y=obj.Y; 168 X = obj.X;
164 xi=obj.xi; eta=obj.eta; 169 Y = obj.Y;
165 170 xi = obj.xi;
166 171 eta = obj.eta;
172
167 switch boundary 173 switch boundary
168 case {'w','W','west'} 174 case {'w','W','west'}
169 e_=obj.e_w; 175 e_ = obj.e_w;
170 mat=obj.Ahat; 176 mat = obj.Ahat;
171 boundPos='l'; 177 boundPos = 'l';
172 Hi=obj.Hxii; 178 Hi = obj.Hxii;
173 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w)); 179 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w));
174 side=max(length(eta)); 180 side = max(length(eta));
175 case {'e','E','east'} 181 case {'e','E','east'}
176 e_=obj.e_e; 182 e_ = obj.e_e;
177 mat=obj.Ahat; 183 mat = obj.Ahat;
178 boundPos='r'; 184 boundPos = 'r';
179 Hi=obj.Hxii; 185 Hi = obj.Hxii;
180 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e)); 186 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e));
181 side=max(length(eta)); 187 side = max(length(eta));
182 case {'s','S','south'} 188 case {'s','S','south'}
183 e_=obj.e_s; 189 e_ = obj.e_s;
184 mat=obj.Bhat; 190 mat = obj.Bhat;
185 boundPos='l'; 191 boundPos = 'l';
186 Hi=obj.Hetai; 192 Hi = obj.Hetai;
187 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s)); 193 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s));
188 side=max(length(xi)); 194 side = max(length(xi));
189 case {'n','N','north'} 195 case {'n','N','north'}
190 e_=obj.e_n; 196 e_ = obj.e_n;
191 mat=obj.Bhat; 197 mat = obj.Bhat;
192 boundPos='r'; 198 boundPos = 'r';
193 Hi=obj.Hetai; 199 Hi = obj.Hetai;
194 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n)); 200 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n));
195 side=max(length(xi)); 201 side = max(length(xi));
196 end 202 end
197 203
198 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 204 pos = signVec(1);
199 205 zeroval = signVec(2);
206 neg = signVec(3);
207
200 switch boundPos 208 switch boundPos
201 case {'l'} 209 case {'l'}
202 tau=sparse(obj.n*side,pos); 210 tau = sparse(obj.n*side,pos);
203 Vi_plus=Vi(1:pos,:); 211 Vi_plus = Vi(1:pos,:);
204 tau(1:pos,:)=-abs(D(1:pos,1:pos)); 212 tau(1:pos,:) = -abs(D(1:pos,1:pos));
205 closure=Hi*e_*V*tau*Vi_plus*e_'; 213 closure = Hi*e_*V*tau*Vi_plus*e_';
206 penalty=-Hi*e_*V*tau*Vi_plus; 214 penalty = -Hi*e_*V*tau*Vi_plus;
207 case {'r'} 215 case {'r'}
208 tau=sparse(obj.n*side,neg); 216 tau = sparse(obj.n*side,neg);
209 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 217 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
210 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:); 218 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
211 closure=Hi*e_*V*tau*Vi_minus*e_'; 219 closure = Hi*e_*V*tau*Vi_minus*e_';
212 penalty=-Hi*e_*V*tau*Vi_minus; 220 penalty = -Hi*e_*V*tau*Vi_minus;
213 end 221 end
214 end 222 end
215 223
216 224
217 function [closure,penalty]=boundary_condition_general(obj,boundary,L) 225 % General boundary condition in the form Lu=g(x)
218 params=obj.params; 226 function [closure,penalty] = boundary_condition_general(obj,boundary,L)
219 X=obj.X; Y=obj.Y; 227 params = obj.params;
220 xi=obj.xi; eta=obj.eta; 228 X = obj.X;
221 229 Y = obj.Y;
230 xi = obj.xi;
231 eta = obj.eta;
232
222 switch boundary 233 switch boundary
223 case {'w','W','west'} 234 case {'w','W','west'}
224 e_=obj.e_w; 235 e_ = obj.e_w;
225 mat=obj.Ahat; 236 mat = obj.Ahat;
226 boundPos='l'; 237 boundPos = 'l';
227 Hi=obj.Hxii; 238 Hi = obj.Hxii;
228 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w)); 239 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w));
229 240
230 Ji_vec=diag(obj.Ji); 241 Ji_vec = diag(obj.Ji);
231 Ji=diag(Ji_vec(obj.index_w)); 242 Ji = diag(Ji_vec(obj.index_w));
232 xi_x=Ji*obj.Y_eta(obj.index_w); 243 xi_x = Ji*obj.Y_eta(obj.index_w);
233 xi_y=-Ji*obj.X_eta(obj.index_w); 244 xi_y = -Ji*obj.X_eta(obj.index_w);
234 L=obj.evaluateCoefficientMatrix(L,xi_x,xi_y,[],[]); 245 L = obj.evaluateCoefficientMatrix(L,xi_x,xi_y,[],[]);
235 side=max(length(eta)); 246 side = max(length(eta));
236 case {'e','E','east'} 247 case {'e','E','east'}
237 e_=obj.e_e; 248 e_ = obj.e_e;
238 mat=obj.Ahat; 249 mat = obj.Ahat;
239 boundPos='r'; 250 boundPos = 'r';
240 Hi=obj.Hxii; 251 Hi = obj.Hxii;
241 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e)); 252 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e));
242 253
243 Ji_vec=diag(obj.Ji); 254 Ji_vec = diag(obj.Ji);
244 Ji=diag(Ji_vec(obj.index_e)); 255 Ji = diag(Ji_vec(obj.index_e));
245 xi_x=Ji*obj.Y_eta(obj.index_e); 256 xi_x = Ji*obj.Y_eta(obj.index_e);
246 xi_y=-Ji*obj.X_eta(obj.index_e); 257 xi_y = -Ji*obj.X_eta(obj.index_e);
247 L=obj.evaluateCoefficientMatrix(L,-xi_x,-xi_y,[],[]); 258 L = obj.evaluateCoefficientMatrix(L,-xi_x,-xi_y,[],[]);
248 side=max(length(eta)); 259 side = max(length(eta));
249 case {'s','S','south'} 260 case {'s','S','south'}
250 e_=obj.e_s; 261 e_ = obj.e_s;
251 mat=obj.Bhat; 262 mat = obj.Bhat;
252 boundPos='l'; 263 boundPos = 'l';
253 Hi=obj.Hetai; 264 Hi = obj.Hetai;
254 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s)); 265 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s));
255 266
256 Ji_vec=diag(obj.Ji); 267 Ji_vec = diag(obj.Ji);
257 Ji=diag(Ji_vec(obj.index_s)); 268 Ji = diag(Ji_vec(obj.index_s));
258 eta_x=Ji*obj.Y_xi(obj.index_s); 269 eta_x = Ji*obj.Y_xi(obj.index_s);
259 eta_y=-Ji*obj.X_xi(obj.index_s); 270 eta_y = -Ji*obj.X_xi(obj.index_s);
260 L=obj.evaluateCoefficientMatrix(L,eta_x,eta_y,[],[]); 271 L = obj.evaluateCoefficientMatrix(L,eta_x,eta_y,[],[]);
261 side=max(length(xi)); 272 side = max(length(xi));
262 case {'n','N','north'} 273 case {'n','N','north'}
263 e_=obj.e_n; 274 e_ = obj.e_n;
264 275 mat = obj.Bhat;
265 mat=obj.Bhat; 276 boundPos = 'r';
266 boundPos='r'; 277 Hi = obj.Hetai;
267 Hi=obj.Hetai; 278 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n));
268 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n)); 279
269 280 Ji_vec = diag(obj.Ji);
270 Ji_vec=diag(obj.Ji); 281 Ji = diag(Ji_vec(obj.index_n));
271 Ji=diag(Ji_vec(obj.index_n)); 282 eta_x = Ji*obj.Y_xi(obj.index_n);
272 eta_x=Ji*obj.Y_xi(obj.index_n); 283 eta_y = -Ji*obj.X_xi(obj.index_n);
273 eta_y=-Ji*obj.X_xi(obj.index_n); 284 L = obj.evaluateCoefficientMatrix(L,-eta_x,-eta_y,[],[]);
274 L=obj.evaluateCoefficientMatrix(L,-eta_x,-eta_y,[],[]); 285 side = max(length(xi));
275 286 end
276 side=max(length(xi)); 287
277 288 pos = signVec(1);
278 end 289 zeroval = signVec(2);
279 290 neg = signVec(3);
280 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 291
281
282 switch boundPos 292 switch boundPos
283 case {'l'} 293 case {'l'}
284 tau=sparse(obj.n*side,pos); 294 tau = sparse(obj.n*side,pos);
285 Vi_plus=Vi(1:pos,:); 295 Vi_plus = Vi(1:pos,:);
286 Vi_minus=Vi(pos+1:obj.n*side,:); 296 Vi_minus = Vi(pos+1:obj.n*side,:);
287 V_plus=V(:,1:pos); 297 V_plus = V(:,1:pos);
288 V_minus=V(:,(pos)+1:obj.n*side); 298 V_minus = V(:,(pos)+1:obj.n*side);
289 299
290 tau(1:pos,:)=-abs(D(1:pos,1:pos)); 300 tau(1:pos,:) = -abs(D(1:pos,1:pos));
291 R=-inv(L*V_plus)*(L*V_minus); 301 R = -inv(L*V_plus)*(L*V_minus);
292 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; 302 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
293 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L; 303 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L;
294 case {'r'} 304 case {'r'}
295 tau=sparse(obj.n*side,neg); 305 tau = sparse(obj.n*side,neg);
296 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 306 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
297 Vi_plus=Vi(1:pos,:); 307 Vi_plus = Vi(1:pos,:);
298 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:); 308 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
299 309
300 V_plus=V(:,1:pos); 310 V_plus = V(:,1:pos);
301 V_minus=V(:,(pos+zeroval)+1:obj.n*side); 311 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
302 R=-inv(L*V_minus)*(L*V_plus); 312 R = -inv(L*V_minus)*(L*V_plus);
303 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 313 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
304 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L; 314 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
305 end 315 end
306 end 316 end
307 317
308 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,x_,y_) 318 % Function that diagonalizes a symbolic matrix A as A=V*D*Vi
309 params=obj.params; 319 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign
310 syms xs ys 320 % [d+ ]
311 if(sum(abs(x_))~=0) 321 % D = [ d0 ]
322 % [ d-]
323 % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D
324 function [V,Vi, D,signVec] = matrixDiag(obj,mat,x,y,x_,y_)
325 params = obj.params;
326 syms xs ys
327 if(sum(abs(x_)) ~= 0)
312 syms xs_ 328 syms xs_
313 else 329 else
314 xs_=0; 330 xs_ = 0;
315 end 331 end
316 332
317 if(sum(abs(y_))~=0) 333 if(sum(abs(y_))~= 0)
318 syms ys_; 334 syms ys_;
319 else 335 else
320 ys_=0; 336 ys_ = 0;
321 end 337 end
322 338
323 [V, D]=eig(mat(params,xs,ys,xs_,ys_)); 339 [V, D] = eig(mat(params,xs,ys,xs_,ys_));
324 Vi=inv(V); 340 Vi = inv(V);
325 syms xs ys xs_ ys_ 341 syms xs ys xs_ ys_
326 342
327 xs=x; 343 xs = x;
328 ys=y; 344 ys = y;
329 xs_=x_; 345 xs_ = x_;
330 ys_=y_; 346 ys_ = y_;
331 347
332 side=max(length(x),length(y)); 348 side = max(length(x),length(y));
333 Dret=zeros(obj.n,side*obj.n); 349 Dret = zeros(obj.n,side*obj.n);
334 Vret=zeros(obj.n,side*obj.n); 350 Vret = zeros(obj.n,side*obj.n);
335 Viret=zeros(obj.n,side*obj.n); 351 Viret = zeros(obj.n,side*obj.n);
336 for ii=1:obj.n 352 for ii = 1:obj.n
337 for jj=1:obj.n 353 for jj = 1:obj.n
338 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii)); 354 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii));
339 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii)); 355 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii));
340 Viret(jj,(ii-1)*side+1:side*ii)=eval(Vi(jj,ii)); 356 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));
341 end 357 end
342 end 358 end
343 359
344 D=sparse(Dret); 360 D = sparse(Dret);
345 V=sparse(Vret); 361 V = sparse(Vret);
346 Vi=sparse(Viret); 362 Vi = sparse(Viret);
347 V=obj.evaluateCoefficientMatrix(V,x,y,x_,y_); 363 V = obj.evaluateCoefficientMatrix(V,x,y,x_,y_);
348 D=obj.evaluateCoefficientMatrix(D,x,y,x_,y_); 364 D = obj.evaluateCoefficientMatrix(D,x,y,x_,y_);
349 Vi=obj.evaluateCoefficientMatrix(Vi,x,y,x_,y_); 365 Vi = obj.evaluateCoefficientMatrix(Vi,x,y,x_,y_);
350 DD=diag(D); 366 DD = diag(D);
351 367
352 poseig=(DD>0); 368 poseig = (DD>0);
353 zeroeig=(DD==0); 369 zeroeig = (DD==0);
354 negeig=(DD<0); 370 negeig = (DD<0);
355 371
356 D=diag([DD(poseig); DD(zeroeig); DD(negeig)]); 372 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]);
357 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)]; 373 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
358 Vi=[Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; 374 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
359 signVec=[sum(poseig),sum(zeroeig),sum(negeig)]; 375 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
360 end 376 end
361 end 377 end
362 end 378 end