Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2dCurve.m @ 1072:6468a5f6ec79 feature/grids/LaplaceSquared

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 12 Feb 2019 17:12:42 +0100
parents 8d73fcdb07a5
children
comparison
equal deleted inserted replaced
1071:92cb03e64ca4 1072:6468a5f6ec79
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 % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu
31 function obj = Hypsyst2dCurve(m, order, A, B, E, params, ti) 31 function obj = Hypsyst2dCurve(m, order, A, B, E, params, ti)
32 default_arg('E', []) 32 default_arg('E', [])
33 xilim = {0 1}; 33 xilim = {0 1};
34 etalim = {0 1}; 34 etalim = {0 1};
35 35
36 if length(m) == 1 36 if length(m) == 1
37 m = [m m]; 37 m = [m m];
38 end 38 end
39 obj.params = params; 39 obj.params = params;
40 obj.A=A; 40 obj.A=A;
41 obj.B=B; 41 obj.B=B;
42 42
43 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);
44 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);
45 obj.E=@(params,x,y,~,~)E(params,x,y); 45 obj.E=@(params,x,y,~,~)E(params,x,y);
46 46
47 m_xi = m(1); 47 m_xi = m(1);
48 m_eta = m(2); 48 m_eta = m(2);
49 m_tot=m_xi*m_eta; 49 m_tot=m_xi*m_eta;
50 50
51 ops_xi = sbp.D2Standard(m_xi,xilim,order); 51 ops_xi = sbp.D2Standard(m_xi,xilim,order);
52 ops_eta = sbp.D2Standard(m_eta,etalim,order); 52 ops_eta = sbp.D2Standard(m_eta,etalim,order);
53 53
54 obj.xi = ops_xi.x; 54 obj.xi = ops_xi.x;
55 obj.eta = ops_eta.x; 55 obj.eta = ops_eta.x;
56 56
57 obj.Xi = kr(obj.xi,ones(m_eta,1)); 57 obj.Xi = kr(obj.xi,ones(m_eta,1));
58 obj.Eta = kr(ones(m_xi,1),obj.eta); 58 obj.Eta = kr(ones(m_xi,1),obj.eta);
59 59
60 obj.n = length(A(obj.params,0,0)); 60 obj.n = length(A(obj.params,0,0));
61 obj.onesN=ones(obj.n); 61 obj.onesN=ones(obj.n);
62 62
63 obj.index_w=1:m_eta; 63 obj.index_w=1:m_eta;
64 obj.index_e=(m_tot-m_e 64 obj.index_e=(m_tot-m_e
65 65
66 metric_termsta+1):m_tot; 66 metric_termsta+1):m_tot;
67 obj.index_s=1:m_eta:(m_tot-m_eta+1); 67 obj.index_s=1:m_eta:(m_tot-m_eta+1);
68 obj.index_n=(m_eta):m_eta:m_tot; 68 obj.index_n=(m_eta):m_eta:m_tot;
69 69
70 I_n = eye(obj.n); 70 I_n = eye(obj.n);
71 I_xi = speye(m_xi); 71 I_xi = speye(m_xi);
72 obj.I_xi = I_xi; 72 obj.I_xi = I_xi;
73 I_eta = speye(m_eta); 73 I_eta = speye(m_eta);
74 obj.I_eta = I_eta; 74 obj.I_eta = I_eta;
75 75
76 D1_xi = kr(I_n, ops_xi.D1, I_eta); 76 D1_xi = kr(I_n, ops_xi.D1, I_eta);
77 obj.Hxii = kr(I_n, ops_xi.HI, I_eta); 77 obj.Hxii = kr(I_n, ops_xi.HI, I_eta);
78 D1_eta = kr(I_n, I_xi, ops_eta.D1); 78 D1_eta = kr(I_n, I_xi, ops_eta.D1);
79 obj.Hetai = kr(I_n, I_xi, ops_eta.HI); 79 obj.Hetai = kr(I_n, I_xi, ops_eta.HI);
80 80
81 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);
82 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);
83 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);
84 obj.e_n = kr(I_n, I_xi, 84 obj.e_n = kr(I_n, I_xi,
85 85
86 metric_termsops_eta.e_r); 86 metric_termsops_eta.e_r);
87 87
88 [X,Y] = ti.map(obj.xi,obj.eta); 88 [X,Y] = ti.map(obj.xi,obj.eta);
89 89
90 [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);
91 [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);
92 92
93 obj.X = reshape(X,m_tot,1); 93 obj.X = reshape(X,m_tot,1);
94 obj.Y = reshape(Y,m_tot,1); 94 obj.Y = reshape(Y,m_tot,1);
95 obj.X_xi = reshape(x_xi,m_tot,1); 95 obj.X_xi = reshape(x_xi,m_tot,1);
96 obj.Y_xi = reshape(y_xi,m_tot,1); 96 obj.Y_xi = reshape(y_xi,m_tot,1);
97 obj.X_eta = reshape(x_eta,m_tot,1); 97 obj.X_eta = reshape(x_eta,m_tot,1);
98 obj.Y_eta = reshape(y_eta,m_tot,1); 98 obj.Y_eta = reshape(y_eta,m_tot,1);
99 99
100 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);
101 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);
102 E_evaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,[],[]); 102 E_evaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,[],[]);
103 103
104 obj.m = m; 104 obj.m = m;
105 obj.h = [ops_xi.h ops_eta.h]; 105 obj.h = [ops_xi.h ops_eta.h];
106 obj.order = order; 106 obj.order = order;
107 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;
108 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));
109 109
110 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;
111 111
112 end 112 end
113 113
114 % 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
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. 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.
116 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w',General boundary conditions'n','s'. 116 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w',General boundary conditions'n','s'.
117 % 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.
118 % 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.
125 [closure,penalty] = boundary_condition_general(obj,boundary,L); 125 [closure,penalty] = boundary_condition_general(obj,boundary,L);
126 otherwise 126 otherwise
127 error('No such boundary condition') 127 error('No such boundary condition')
128 end 128 end
129 end 129 end
130 130
131 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundaryGeneral boundary conditions) 131 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
132 error('An interface function does not exist yet'); 132 error('Not implemented');
133 end 133 end
134 134
135 function N = size(obj) 135 function N = size(obj)
136 N = obj.m; 136 N = obj.m;
137 end 137 end
138 138
139 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_) 139 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_)
140 params = obj.params; 140 params = obj.params;
141 141
142 if isa(mat,'function_handle') 142 if isa(mat,'function_handle')
143 [rows,cols] = size(mat(params,0,0,0,0)); 143 [rows,cols] = size(mat(params,0,0,0,0));
144 x_ = kr(obj.onesN,x_); 144 x_ = kr(obj.onesN,x_);
145 y_ = kr(obj.onesN,y_); 145 y_ = kr(obj.onesN,y_);
146 matVec = mat(params,X',Y',x_',y_'); 146 matVec = mat(params,X',Y',x_',y_');
150 matVec = mat; 150 matVec = mat;
151 [rows,cols] = size(matVec); 151 [rows,cols] = size(matVec);
152 side = max(length(X),length(Y)); 152 side = max(length(X),length(Y));
153 cols = cols/side; 153 cols = cols/side;
154 end 154 end
155 155
156 ret = cell(rows,cols); 156 ret = cell(rows,cols);
157 for ii = 1:rows 157 for ii = 1:rows
158 for jj = 1:cols 158 for jj = 1:cols
159 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));
160 end 160 end
161 end 161 end
162 ret = cell2mat(ret); 162 ret = cell2mat(ret);
163 end 163 end
164 164
165 %Characteristic boundary conditions 165 %Characteristic boundary conditions
166 function [closure, penalty] = boundary_condition_char(obj,boundary) 166 function [closure, penalty] = boundary_condition_char(obj,boundary)
167 params = obj.params; 167 params = obj.params;
168 X = obj.X; 168 X = obj.X;
169 Y = obj.Y; 169 Y = obj.Y;
170 xi = obj.xi; 170 xi = obj.xi;
171 eta = obj.eta; 171 eta = obj.eta;
172 172 e_ = obj.getBoundaryOperator('e', boundary);
173
173 switch boundary 174 switch boundary
174 case {'w','W','west'} 175 case {'w','W','west'}
175 e_ = obj.e_w;
176 mat = obj.Ahat; 176 mat = obj.Ahat;
177 boundPos = 'l'; 177 boundPos = 'l';
178 Hi = obj.Hxii; 178 Hi = obj.Hxii;
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)); 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));
180 side = max(length(eta)); 180 side = max(length(eta));
181 case {'e','E','east'} 181 case {'e','E','east'}
182 e_ = obj.e_e;
183 mat = obj.Ahat; 182 mat = obj.Ahat;
184 boundPos = 'r'; 183 boundPos = 'r';
185 Hi = obj.Hxii; 184 Hi = obj.Hxii;
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)); 185 [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));
187 side = max(length(eta)); 186 side = max(length(eta));
188 case {'s','S','south'} 187 case {'s','S','south'}
189 e_ = obj.e_s;
190 mat = obj.Bhat; 188 mat = obj.Bhat;
191 boundPos = 'l'; 189 boundPos = 'l';
192 Hi = obj.Hetai; 190 Hi = obj.Hetai;
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)); 191 [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));
194 side = max(length(xi)); 192 side = max(length(xi));
195 case {'n','N','north'} 193 case {'n','N','north'}
196 e_ = obj.e_n;
197 mat = obj.Bhat; 194 mat = obj.Bhat;
198 boundPos = 'r'; 195 boundPos = 'r';
199 Hi = obj.Hetai; 196 Hi = obj.Hetai;
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)); 197 [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));
201 side = max(length(xi)); 198 side = max(length(xi));
202 end 199 end
203 200
204 pos = signVec(1); 201 pos = signVec(1);
205 zeroval = signVec(2); 202 zeroval = signVec(2);
206 neg = signVec(3); 203 neg = signVec(3);
207 204
208 switch boundPos 205 switch boundPos
209 case {'l'} 206 case {'l'}
210 tau = sparse(obj.n*side,pos); 207 tau = sparse(obj.n*side,pos);
211 Vi_plus = Vi(1:pos,:); 208 Vi_plus = Vi(1:pos,:);
212 tau(1:pos,:) = -abs(D(1:pos,1:pos)); 209 tau(1:pos,:) = -abs(D(1:pos,1:pos));
216 tau = sparse(obj.n*side,neg); 213 tau = sparse(obj.n*side,neg);
217 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 214 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
218 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); 215 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
219 closure = Hi*e_*V*tau*Vi_minus*e_'; 216 closure = Hi*e_*V*tau*Vi_minus*e_';
220 penalty = -Hi*e_*V*tau*Vi_minus; 217 penalty = -Hi*e_*V*tau*Vi_minus;
221 end 218 end
222 end 219 end
223 220
224 221
225 % General boundary condition in the form Lu=g(x) 222 % General boundary condition in the form Lu=g(x)
226 function [closure,penalty] = boundary_condition_general(obj,boundary,L) 223 function [closure,penalty] = boundary_condition_general(obj,boundary,L)
227 params = obj.params; 224 params = obj.params;
228 X = obj.X; 225 X = obj.X;
229 Y = obj.Y; 226 Y = obj.Y;
230 xi = obj.xi; 227 xi = obj.xi;
231 eta = obj.eta; 228 eta = obj.eta;
232 229
233 switch boundary 230 switch boundary
234 case {'w','W','west'} 231 case {'w','W','west'}
235 e_ = obj.e_w; 232 e_ = obj.e_w;
236 mat = obj.Ahat; 233 mat = obj.Ahat;
237 boundPos = 'l'; 234 boundPos = 'l';
238 Hi = obj.Hxii; 235 Hi = obj.Hxii;
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)); 236 [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));
240 237
241 Ji_vec = diag(obj.Ji); 238 Ji_vec = diag(obj.Ji);
242 Ji = diag(Ji_vec(obj.index_w)); 239 Ji = diag(Ji_vec(obj.index_w));
243 xi_x = Ji*obj.Y_eta(obj.index_w); 240 xi_x = Ji*obj.Y_eta(obj.index_w);
244 xi_y = -Ji*obj.X_eta(obj.index_w); 241 xi_y = -Ji*obj.X_eta(obj.index_w);
245 L = obj.evaluateCoefficientMatrix(L,xi_x,xi_y,[],[]); 242 L = obj.evaluateCoefficientMatrix(L,xi_x,xi_y,[],[]);
248 e_ = obj.e_e; 245 e_ = obj.e_e;
249 mat = obj.Ahat; 246 mat = obj.Ahat;
250 boundPos = 'r'; 247 boundPos = 'r';
251 Hi = obj.Hxii; 248 Hi = obj.Hxii;
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)); 249 [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));
253 250
254 Ji_vec = diag(obj.Ji); 251 Ji_vec = diag(obj.Ji);
255 Ji = diag(Ji_vec(obj.index_e)); 252 Ji = diag(Ji_vec(obj.index_e));
256 xi_x = Ji*obj.Y_eta(obj.index_e); 253 xi_x = Ji*obj.Y_eta(obj.index_e);
257 xi_y = -Ji*obj.X_eta(obj.index_e); 254 xi_y = -Ji*obj.X_eta(obj.index_e);
258 L = obj.evaluateCoefficientMatrix(L,-xi_x,-xi_y,[],[]); 255 L = obj.evaluateCoefficientMatrix(L,-xi_x,-xi_y,[],[]);
261 e_ = obj.e_s; 258 e_ = obj.e_s;
262 mat = obj.Bhat; 259 mat = obj.Bhat;
263 boundPos = 'l'; 260 boundPos = 'l';
264 Hi = obj.Hetai; 261 Hi = obj.Hetai;
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)); 262 [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));
266 263
267 Ji_vec = diag(obj.Ji); 264 Ji_vec = diag(obj.Ji);
268 Ji = diag(Ji_vec(obj.index_s)); 265 Ji = diag(Ji_vec(obj.index_s));
269 eta_x = Ji*obj.Y_xi(obj.index_s); 266 eta_x = Ji*obj.Y_xi(obj.index_s);
270 eta_y = -Ji*obj.X_xi(obj.index_s); 267 eta_y = -Ji*obj.X_xi(obj.index_s);
271 L = obj.evaluateCoefficientMatrix(L,eta_x,eta_y,[],[]); 268 L = obj.evaluateCoefficientMatrix(L,eta_x,eta_y,[],[]);
274 e_ = obj.e_n; 271 e_ = obj.e_n;
275 mat = obj.Bhat; 272 mat = obj.Bhat;
276 boundPos = 'r'; 273 boundPos = 'r';
277 Hi = obj.Hetai; 274 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)); 275 [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 276
280 Ji_vec = diag(obj.Ji); 277 Ji_vec = diag(obj.Ji);
281 Ji = diag(Ji_vec(obj.index_n)); 278 Ji = diag(Ji_vec(obj.index_n));
282 eta_x = Ji*obj.Y_xi(obj.index_n); 279 eta_x = Ji*obj.Y_xi(obj.index_n);
283 eta_y = -Ji*obj.X_xi(obj.index_n); 280 eta_y = -Ji*obj.X_xi(obj.index_n);
284 L = obj.evaluateCoefficientMatrix(L,-eta_x,-eta_y,[],[]); 281 L = obj.evaluateCoefficientMatrix(L,-eta_x,-eta_y,[],[]);
285 side = max(length(xi)); 282 side = max(length(xi));
286 end 283 end
287 284
288 pos = signVec(1); 285 pos = signVec(1);
289 zeroval = signVec(2); 286 zeroval = signVec(2);
290 neg = signVec(3); 287 neg = signVec(3);
291 288
292 switch boundPos 289 switch boundPos
293 case {'l'} 290 case {'l'}
294 tau = sparse(obj.n*side,pos); 291 tau = sparse(obj.n*side,pos);
295 Vi_plus = Vi(1:pos,:); 292 Vi_plus = Vi(1:pos,:);
296 Vi_minus = Vi(pos+1:obj.n*side,:); 293 Vi_minus = Vi(pos+1:obj.n*side,:);
297 V_plus = V(:,1:pos); 294 V_plus = V(:,1:pos);
298 V_minus = V(:,(pos)+1:obj.n*side); 295 V_minus = V(:,(pos)+1:obj.n*side);
299 296
300 tau(1:pos,:) = -abs(D(1:pos,1:pos)); 297 tau(1:pos,:) = -abs(D(1:pos,1:pos));
301 R = -inv(L*V_plus)*(L*V_minus); 298 R = -inv(L*V_plus)*(L*V_minus);
302 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; 299 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
303 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L; 300 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L;
304 case {'r'} 301 case {'r'}
305 tau = sparse(obj.n*side,neg); 302 tau = sparse(obj.n*side,neg);
306 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 303 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
307 Vi_plus = Vi(1:pos,:); 304 Vi_plus = Vi(1:pos,:);
308 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); 305 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
309 306
310 V_plus = V(:,1:pos); 307 V_plus = V(:,1:pos);
311 V_minus = V(:,(pos+zeroval)+1:obj.n*side); 308 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
312 R = -inv(L*V_minus)*(L*V_plus); 309 R = -inv(L*V_minus)*(L*V_plus);
313 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 310 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
314 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; 311 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
315 end 312 end
316 end 313 end
317 314
318 % Function that diagonalizes a symbolic matrix A as A=V*D*Vi 315 % Function that diagonalizes a symbolic matrix A as A=V*D*Vi
319 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign 316 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign
320 % [d+ ] 317 % [d+ ]
321 % D = [ d0 ] 318 % D = [ d0 ]
322 % [ d-] 319 % [ d-]
327 if(sum(abs(x_)) ~= 0) 324 if(sum(abs(x_)) ~= 0)
328 syms xs_ 325 syms xs_
329 else 326 else
330 xs_ = 0; 327 xs_ = 0;
331 end 328 end
332 329
333 if(sum(abs(y_))~= 0) 330 if(sum(abs(y_))~= 0)
334 syms ys_; 331 syms ys_;
335 else 332 else
336 ys_ = 0; 333 ys_ = 0;
337 end 334 end
338 335
339 [V, D] = eig(mat(params,xs,ys,xs_,ys_)); 336 [V, D] = eig(mat(params,xs,ys,xs_,ys_));
340 Vi = inv(V); 337 Vi = inv(V);
341 syms xs ys xs_ ys_ 338 syms xs ys xs_ ys_
342 339
343 xs = x; 340 xs = x;
344 ys = y; 341 ys = y;
345 xs_ = x_; 342 xs_ = x_;
346 ys_ = y_; 343 ys_ = y_;
347 344
348 side = max(length(x),length(y)); 345 side = max(length(x),length(y));
349 Dret = zeros(obj.n,side*obj.n); 346 Dret = zeros(obj.n,side*obj.n);
350 Vret = zeros(obj.n,side*obj.n); 347 Vret = zeros(obj.n,side*obj.n);
351 Viret = zeros(obj.n,side*obj.n); 348 Viret = zeros(obj.n,side*obj.n);
352 for ii = 1:obj.n 349 for ii = 1:obj.n
354 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); 351 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii));
355 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii)); 352 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii));
356 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); 353 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));
357 end 354 end
358 end 355 end
359 356
360 D = sparse(Dret); 357 D = sparse(Dret);
361 V = sparse(Vret); 358 V = sparse(Vret);
362 Vi = sparse(Viret); 359 Vi = sparse(Viret);
363 V = obj.evaluateCoefficientMatrix(V,x,y,x_,y_); 360 V = obj.evaluateCoefficientMatrix(V,x,y,x_,y_);
364 D = obj.evaluateCoefficientMatrix(D,x,y,x_,y_); 361 D = obj.evaluateCoefficientMatrix(D,x,y,x_,y_);
365 Vi = obj.evaluateCoefficientMatrix(Vi,x,y,x_,y_); 362 Vi = obj.evaluateCoefficientMatrix(Vi,x,y,x_,y_);
366 DD = diag(D); 363 DD = diag(D);
367 364
368 poseig = (DD>0); 365 poseig = (DD>0);
369 zeroeig = (DD==0); 366 zeroeig = (DD==0);
370 negeig = (DD<0); 367 negeig = (DD<0);
371 368
372 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); 369 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]);
373 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; 370 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
374 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; 371 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
375 signVec = [sum(poseig),sum(zeroeig),sum(negeig)]; 372 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
376 end 373 end
374
375 % Returns the boundary operator op for the boundary specified by the string boundary.
376 % op -- string or a cell array of strings
377 % boundary -- string
378 function varargout = getBoundaryOperator(obj, op, boundary)
379 assertIsMember(boundary, {'w', 'e', 's', 'n'})
380
381 if ~iscell(op)
382 op = {op};
383 end
384
385 for i = 1:numel(op)
386 switch op{i}
387 case 'e'
388 switch boundary
389 case 'w'
390 e = obj.e_w;
391 case 'e'
392 e = obj.e_e;
393 case 's'
394 e = obj.e_s;
395 case 'n'
396 e = obj.e_n;
397 end
398 varargout{i} = e;
399 end
400 end
401 end
402
403 % Returns square boundary quadrature matrix, of dimension
404 % corresponding to the number of boundary points
405 %
406 % boundary -- string
407 function H_b = getBoundaryQuadrature(obj, boundary)
408 assertIsMember(boundary, {'w', 'e', 's', 'n'})
409
410 e = obj.getBoundaryOperator('e', boundary);
411
412 switch boundary
413 case 'w'
414 H_b = inv(e'*obj.Hetai*e);
415 case 'e'
416 H_b = inv(e'*obj.Hetai*e);
417 case 's'
418 H_b = inv(e'*obj.Hxii*e);
419 case 'n'
420 H_b = inv(e'*obj.Hxii*e);
421 end
422 end
423
424
377 end 425 end
378 end 426 end