Mercurial > repos > public > sbplib

comparison +scheme/Wave2dCurve.m @ 820:501750fbbfdb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/grids
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 07 Sep 2018 14:40:58 +0200
parents 4e266dfe9edc
children 459eeb99130f
comparison
equal deleted inserted replaced
819:fdf0ef9150f4 820:501750fbbfdb
1 classdef Wave2dCurve < scheme.Scheme 1 classdef Wave2dCurve < 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 h % Grid spacing 4 h % Grid spacing
5 u,v % Grid 5
6 x,y % Values of x and y for each grid point 6 grid
7 X,Y % Grid point locations as matrices 7
8 order % Order accuracy for the approximation 8 order % Order accuracy for the approximation
9 9
10 D % non-stabalized scheme operator 10 D % non-stabalized scheme operator
11 M % Derivative norm 11 M % Derivative norm
12 c 12 c
24 du_e, dv_e 24 du_e, dv_e
25 du_s, dv_s 25 du_s, dv_s
26 du_n, dv_n 26 du_n, dv_n
27 gamm_u, gamm_v 27 gamm_u, gamm_v
28 lambda 28 lambda
29
30 Dx, Dy % Physical derivatives
31
32 x_u
33 x_v
34 y_u
35 y_v
29 end 36 end
30 37
31 methods 38 methods
32 function obj = Wave2dCurve(m,ti,order,c,opSet) 39 function obj = Wave2dCurve(g ,order, c, opSet)
33 default_arg('opSet',@sbp.D2Variable); 40 default_arg('opSet',@sbp.D2Variable);
34 default_arg('c', 1); 41 default_arg('c', 1);
35 42
36 if length(m) == 1 43 warning('Use LaplaceCruveilinear instead')
37 m = [m m]; 44
38 end 45 assert(isa(g, 'grid.Curvilinear'))
39 46
47 m = g.size();
40 m_u = m(1); 48 m_u = m(1);
41 m_v = m(2); 49 m_v = m(2);
42 m_tot = m_u*m_v; 50 m_tot = g.N();
43 51
44 [u, h_u] = util.get_grid(0, 1, m_u); 52 h = g.scaling();
45 [v, h_v] = util.get_grid(0, 1, m_v); 53 h_u = h(1);
46 54 h_v = h(2);
47 55
48 % Operators 56 % Operators
49 ops_u = opSet(m_u, {0, 1}, order); 57 ops_u = opSet(m_u, {0, 1}, order);
50 ops_v = opSet(m_v, {0, 1}, order); 58 ops_v = opSet(m_v, {0, 1}, order);
51 59
52 I_u = speye(m_u); 60 I_u = speye(m_u);
53 I_v = speye(m_v); 61 I_v = speye(m_v);
54 62
55 D1_u = sparse(ops_u.D1); 63 D1_u = ops_u.D1;
56 D2_u = ops_u.D2; 64 D2_u = ops_u.D2;
57 H_u = sparse(ops_u.H); 65 H_u = ops_u.H;
58 Hi_u = sparse(ops_u.HI); 66 Hi_u = ops_u.HI;
59 % M_u = sparse(ops_u.M); 67 e_l_u = ops_u.e_l;
60 e_l_u = sparse(ops_u.e_l); 68 e_r_u = ops_u.e_r;
61 e_r_u = sparse(ops_u.e_r); 69 d1_l_u = ops_u.d1_l;
62 d1_l_u = sparse(ops_u.d1_l); 70 d1_r_u = ops_u.d1_r;
63 d1_r_u = sparse(ops_u.d1_r); 71
64 72 D1_v = ops_v.D1;
65 D1_v = sparse(ops_v.D1);
66 D2_v = ops_v.D2; 73 D2_v = ops_v.D2;
67 H_v = sparse(ops_v.H); 74 H_v = ops_v.H;
68 Hi_v = sparse(ops_v.HI); 75 Hi_v = ops_v.HI;
69 % M_v = sparse(ops_v.M); 76 e_l_v = ops_v.e_l;
70 e_l_v = sparse(ops_v.e_l); 77 e_r_v = ops_v.e_r;
71 e_r_v = sparse(ops_v.e_r); 78 d1_l_v = ops_v.d1_l;
72 d1_l_v = sparse(ops_v.d1_l); 79 d1_r_v = ops_v.d1_r;
73 d1_r_v = sparse(ops_v.d1_r); 80
81 Du = kr(D1_u,I_v);
82 Dv = kr(I_u,D1_v);
74 83
75 % Metric derivatives 84 % Metric derivatives
76 [X,Y] = ti.map(u,v); 85 coords = g.points();
77 86 x = coords(:,1);
78 [x_u,x_v] = gridDerivatives(X,D1_u,D1_v); 87 y = coords(:,2);
79 [y_u,y_v] = gridDerivatives(Y,D1_u,D1_v); 88
80 89 x_u = Du*x;
81 90 x_v = Dv*x;
91 y_u = Du*y;
92 y_v = Dv*y;
82 93
83 J = x_u.*y_v - x_v.*y_u; 94 J = x_u.*y_v - x_v.*y_u;
84 a11 = 1./J .* (x_v.^2 + y_v.^2); %% GÖR SOM MATRISER 95 a11 = 1./J .* (x_v.^2 + y_v.^2);
85 a12 = -1./J .* (x_u.*x_v + y_u.*y_v); 96 a12 = -1./J .* (x_u.*x_v + y_u.*y_v);
86 a22 = 1./J .* (x_u.^2 + y_u.^2); 97 a22 = 1./J .* (x_u.^2 + y_u.^2);
87 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2)); 98 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2));
88 99
89 dof_order = reshape(1:m_u*m_v,m_v,m_u); 100 % Assemble full operators
101 L_12 = spdiags(a12, 0, m_tot, m_tot);
102 Duv = Du*L_12*Dv;
103 Dvu = Dv*L_12*Du;
90 104
91 Duu = sparse(m_tot); 105 Duu = sparse(m_tot);
92 Dvv = sparse(m_tot); 106 Dvv = sparse(m_tot);
107 ind = grid.funcToMatrix(g, 1:m_tot);
93 108
94 for i = 1:m_v 109 for i = 1:m_v
95 D = D2_u(a11(i,:)); 110 D = D2_u(a11(ind(:,i)));
96 p = dof_order(i,:); 111 p = ind(:,i);
97 Duu(p,p) = D; 112 Duu(p,p) = D;
98 end 113 end
99 114
100 for i = 1:m_u 115 for i = 1:m_u
101 D = D2_v(a22(:,i)); 116 D = D2_v(a22(ind(i,:)));
102 p = dof_order(:,i); 117 p = ind(i,:);
103 Dvv(p,p) = D; 118 Dvv(p,p) = D;
104 end 119 end
105
106 L_12 = spdiags(a12(:),0,m_tot,m_tot);
107 Du = kr(D1_u,I_v);
108 Dv = kr(I_u,D1_v);
109
110 Duv = Du*L_12*Dv;
111 Dvu = Dv*L_12*Du;
112
113
114 120
115 obj.H = kr(H_u,H_v); 121 obj.H = kr(H_u,H_v);
116 obj.Hi = kr(Hi_u,Hi_v); 122 obj.Hi = kr(Hi_u,Hi_v);
117 obj.Hu = kr(H_u,I_v); 123 obj.Hu = kr(H_u,I_v);
118 obj.Hv = kr(I_u,H_v); 124 obj.Hv = kr(I_u,H_v);
119 obj.Hiu = kr(Hi_u,I_v); 125 obj.Hiu = kr(Hi_u,I_v);
120 obj.Hiv = kr(I_u,Hi_v); 126 obj.Hiv = kr(I_u,Hi_v);
121 127
122 % obj.M = kr(M_u,H_v)+kr(H_u,M_v);
123 obj.e_w = kr(e_l_u,I_v); 128 obj.e_w = kr(e_l_u,I_v);
124 obj.e_e = kr(e_r_u,I_v); 129 obj.e_e = kr(e_r_u,I_v);
125 obj.e_s = kr(I_u,e_l_v); 130 obj.e_s = kr(I_u,e_l_v);
126 obj.e_n = kr(I_u,e_r_v); 131 obj.e_n = kr(I_u,e_r_v);
127 obj.du_w = kr(d1_l_u,I_v); 132 obj.du_w = kr(d1_l_u,I_v);
131 obj.du_s = (obj.e_s'*Du)'; 136 obj.du_s = (obj.e_s'*Du)';
132 obj.dv_s = kr(I_u,d1_l_v); 137 obj.dv_s = kr(I_u,d1_l_v);
133 obj.du_n = (obj.e_n'*Du)'; 138 obj.du_n = (obj.e_n'*Du)';
134 obj.dv_n = kr(I_u,d1_r_v); 139 obj.dv_n = kr(I_u,d1_r_v);
135 140
141 obj.x_u = x_u;
142 obj.x_v = x_v;
143 obj.y_u = y_u;
144 obj.y_v = y_v;
145
136 obj.m = m; 146 obj.m = m;
137 obj.h = [h_u h_v]; 147 obj.h = [h_u h_v];
138 obj.order = order; 148 obj.order = order;
139 149 obj.grid = g;
140 150
141 obj.c = c; 151 obj.c = c;
142 obj.J = spdiags(J(:),0,m_tot,m_tot); 152 obj.J = spdiags(J, 0, m_tot, m_tot);
143 obj.Ji = spdiags(1./J(:),0,m_tot,m_tot); 153 obj.Ji = spdiags(1./J, 0, m_tot, m_tot);
144 obj.a11 = a11; 154 obj.a11 = a11;
145 obj.a12 = a12; 155 obj.a12 = a12;
146 obj.a22 = a22; 156 obj.a22 = a22;
147 obj.D = obj.Ji*c^2*(Duu + Duv + Dvu + Dvv); 157 obj.D = obj.Ji*c^2*(Duu + Duv + Dvu + Dvv);
148 obj.u = u;
149 obj.v = v;
150 obj.X = X;
151 obj.Y = Y;
152 obj.x = X(:);
153 obj.y = Y(:);
154 obj.lambda = lambda; 158 obj.lambda = lambda;
155 159
156 obj.gamm_u = h_u*ops_u.borrowing.M.S; 160 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
157 obj.gamm_v = h_v*ops_v.borrowing.M.S; 161 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
162
163 obj.gamm_u = h_u*ops_u.borrowing.M.d1;
164 obj.gamm_v = h_v*ops_v.borrowing.M.d1;
158 end 165 end
159 166
160 167
161 % Closure functions return the opertors applied to the own doamin to close the boundary 168 % Closure functions return the opertors applied to the own doamin to close the boundary
162 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 169 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
163 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 170 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
164 % type is a string specifying the type of boundary condition if there are several. 171 % type is a string specifying the type of boundary condition if there are several.
165 % data is a function returning the data that should be applied at the boundary. 172 % data is a function returning the data that should be applied at the boundary.
166 % neighbour_scheme is an instance of Scheme that should be interfaced to. 173 % neighbour_scheme is an instance of Scheme that should be interfaced to.
167 % neighbour_boundary is a string specifying which boundary to interface to. 174 % neighbour_boundary is a string specifying which boundary to interface to.
168 function [closure, penalty] = boundary_condition(obj,boundary,type,data) 175 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
169 default_arg('type','neumann'); 176 default_arg('type','neumann');
170 default_arg('data',0); 177 default_arg('parameter', []);
171 178
172 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv] = obj.get_boundary_ops(boundary); 179 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv , ~, ~, ~, scale_factor] = obj.get_boundary_ops(boundary);
173
174 switch type 180 switch type
175 % Dirichlet boundary condition 181 % Dirichlet boundary condition
176 case {'D','d','dirichlet'} 182 case {'D','d','dirichlet'}
177 % v denotes the solution in the neighbour domain 183 % v denotes the solution in the neighbour domain
178 tuning = 1.2; 184 tuning = 1.2;
188 194
189 b1 = gamm*u.lambda./u.a11.^2; 195 b1 = gamm*u.lambda./u.a11.^2;
190 b2 = gamm*u.lambda./u.a22.^2; 196 b2 = gamm*u.lambda./u.a22.^2;
191 197
192 tau = -1./b1 - 1./b2; 198 tau = -1./b1 - 1./b2;
193 tau = tuning * spdiag(tau(:)); 199 tau = tuning * spdiag(tau);
194 sig1 = 1/2; 200 sig1 = 1;
195 201
196 penalty_parameter_1 = halfnorm_inv_n*(tau + sig1*halfnorm_inv_t*F*e'*halfnorm_t)*e; 202 penalty_parameter_1 = halfnorm_inv_n*(tau + sig1*halfnorm_inv_t*F*e'*halfnorm_t)*e;
197 203
198 closure = obj.Ji*obj.c^2 * penalty_parameter_1*e'; 204 closure = obj.Ji*obj.c^2 * penalty_parameter_1*e';
199 pp = -obj.Ji*obj.c^2 * penalty_parameter_1; 205 penalty = -obj.Ji*obj.c^2 * penalty_parameter_1;
200 switch class(data)
201 case 'double'
202 penalty = pp*data;
203 case 'function_handle'
204 penalty = @(t)pp*data(t);
205 otherwise
206 error('Weird data argument!')
207 end
208 206
209 207
210 % Neumann boundary condition 208 % Neumann boundary condition
211 case {'N','n','neumann'} 209 case {'N','n','neumann'}
212 c = obj.c; 210 c = obj.c;
213 211
214
215 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); 212 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n));
216 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); 213 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t));
217 d = (a_n * d_n' + a_t*d_t')'; 214 d = (a_n * d_n' + a_t*d_t')';
218 215
219 tau1 = -s; 216 tau1 = -s;
220 tau2 = 0; 217 tau2 = 0;
221 tau = c.^2 * obj.Ji*(tau1*e + tau2*d); 218 tau = c.^2 * obj.Ji*(tau1*e + tau2*d);
222 219
223 closure = halfnorm_inv*tau*d'; 220 closure = halfnorm_inv*tau*d';
224 221 penalty = -halfnorm_inv*tau;
225 pp = halfnorm_inv*tau; 222
226 switch class(data) 223 % Characteristic boundary condition
227 case 'double' 224 case {'characteristic', 'char', 'c'}
228 penalty = pp*data; 225 default_arg('parameter', 1);
229 case 'function_handle' 226 beta = parameter;
230 penalty = @(t)pp*data(t); 227 c = obj.c;
231 otherwise 228
232 error('Weird data argument!') 229 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n));
233 end 230 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t));
231 d = s*(a_n * d_n' + a_t*d_t')'; % outward facing normal derivative
232
233 tau = -c.^2 * 1/beta*obj.Ji*e;
234
235 warning('is this right?! /c?')
236 closure{1} = halfnorm_inv*tau/c*spdiag(scale_factor)*e';
237 closure{2} = halfnorm_inv*tau*beta*d';
238 penalty = -halfnorm_inv*tau;
234 239
235 % Unknown, boundary condition 240 % Unknown, boundary condition
236 otherwise 241 otherwise
237 error('No such boundary condition: type = %s',type); 242 error('No such boundary condition: type = %s',type);
238 end 243 end
261 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2; 266 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
262 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; 267 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
263 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2; 268 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
264 269
265 tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); 270 tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v);
266 tau = tuning * spdiag(tau(:)); 271 tau = tuning * spdiag(tau);
267 sig1 = 1/2; 272 sig1 = 1/2;
268 sig2 = -1/2; 273 sig2 = -1/2;
269 274
270 penalty_parameter_1 = halfnorm_inv_u_n*(e_u*tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u); 275 penalty_parameter_1 = halfnorm_inv_u_n*(e_u*tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u);
271 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u; 276 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u;
277 282
278 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 283 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
279 % The right boundary is considered the positive boundary 284 % The right boundary is considered the positive boundary
280 % 285 %
281 % I -- the indecies of the boundary points in the grid matrix 286 % I -- the indecies of the boundary points in the grid matrix
282 function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I] = get_boundary_ops(obj,boundary) 287 function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor] = get_boundary_ops(obj, boundary)
283 288
284 gridMatrix = zeros(obj.m(2),obj.m(1)); 289 % gridMatrix = zeros(obj.m(2),obj.m(1));
285 gridMatrix(:) = 1:numel(gridMatrix); 290 % gridMatrix(:) = 1:numel(gridMatrix);
291
292 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
286 293
287 switch boundary 294 switch boundary
288 case 'w' 295 case 'w'
289 e = obj.e_w; 296 e = obj.e_w;
290 d_n = obj.du_w; 297 d_n = obj.du_w;
291 d_t = obj.dv_w; 298 d_t = obj.dv_w;
292 s = -1; 299 s = -1;
293 300
294 I = gridMatrix(:,1); 301 I = ind(1,:);
295 coeff_n = obj.a11(I); 302 coeff_n = obj.a11(I);
296 coeff_t = obj.a12(I); 303 coeff_t = obj.a12(I);
304 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
297 case 'e' 305 case 'e'
298 e = obj.e_e; 306 e = obj.e_e;
299 d_n = obj.du_e; 307 d_n = obj.du_e;
300 d_t = obj.dv_e; 308 d_t = obj.dv_e;
301 s = 1; 309 s = 1;
302 310
303 I = gridMatrix(:,end); 311 I = ind(end,:);
304 coeff_n = obj.a11(I); 312 coeff_n = obj.a11(I);
305 coeff_t = obj.a12(I); 313 coeff_t = obj.a12(I);
314 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
306 case 's' 315 case 's'
307 e = obj.e_s; 316 e = obj.e_s;
308 d_n = obj.dv_s; 317 d_n = obj.dv_s;
309 d_t = obj.du_s; 318 d_t = obj.du_s;
310 s = -1; 319 s = -1;
311 320
312 I = gridMatrix(1,:)'; 321 I = ind(:,1)';
313 coeff_n = obj.a22(I); 322 coeff_n = obj.a22(I);
314 coeff_t = obj.a12(I); 323 coeff_t = obj.a12(I);
324 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
315 case 'n' 325 case 'n'
316 e = obj.e_n; 326 e = obj.e_n;
317 d_n = obj.dv_n; 327 d_n = obj.dv_n;
318 d_t = obj.du_n; 328 d_t = obj.du_n;
319 s = 1; 329 s = 1;
320 330
321 I = gridMatrix(end,:)'; 331 I = ind(:,end)';
322 coeff_n = obj.a22(I); 332 coeff_n = obj.a22(I);
323 coeff_t = obj.a12(I); 333 coeff_t = obj.a12(I);
334 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
324 otherwise 335 otherwise
325 error('No such boundary: boundary = %s',boundary); 336 error('No such boundary: boundary = %s',boundary);
326 end 337 end
327 338
328 switch boundary 339 switch boundary
341 352
342 function N = size(obj) 353 function N = size(obj)
343 N = prod(obj.m); 354 N = prod(obj.m);
344 end 355 end
345 356
346 end 357
347
348 methods(Static)
349 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
350 % and bound_v of scheme schm_v.
351 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
352 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
353 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
354 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
355 end
356 end 358 end
357 end 359 end