Mercurial > repos > public > sbplib

comparison +scheme/Wave2dCurve.m @ 425:e56dbd9e4196 feature/grids

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