Mercurial > repos > public > sbplib

comparison +scheme/Wave2dCurve.m @ 337:e070ebd94d9d feature/beams

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