Mercurial > repos > public > sbplib
annotate +scheme/Wave2dCurve.m @ 87:0a29a60e0b21
In Curve: Rearranged for speed. arc_length_fun is now a property of Curve. If it is not supplied, it is computed via the derivative and spline fitting. Switching to the arc length parameterization is much faster now. The new stuff can be tested with testArcLength.m (which should be deleted after that).
| author | Martin Almquist <martin.almquist@it.uu.se> |
|---|---|
| date | Sun, 29 Nov 2015 22:23:09 +0100 |
| parents | a8ed986fcf57 |
| children | df642750e8f7 |
| rev | line source |
|---|---|
| 0 | 1 classdef Wave2dCurve < scheme.Scheme |
| 2 properties | |
| 3 m % Number of points in each direction, possibly a vector | |
| 4 h % Grid spacing | |
| 5 u,v % Grid | |
| 6 x,y % Values of x and y for each grid point | |
| 7 X,Y % Grid point locations as matrices | |
| 8 order % Order accuracy for the approximation | |
| 9 | |
| 10 D % non-stabalized scheme operator | |
| 11 M % Derivative norm | |
| 12 c | |
| 13 J, Ji | |
| 14 a11, a12, a22 | |
| 15 | |
| 16 H % Discrete norm | |
| 17 Hi | |
| 18 H_u, H_v % Norms in the x and y directions | |
| 19 Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. | |
| 20 Hi_u, Hi_v | |
| 21 Hiu, Hiv | |
| 22 e_w, e_e, e_s, e_n | |
| 23 du_w, dv_w | |
| 24 du_e, dv_e | |
| 25 du_s, dv_s | |
| 26 du_n, dv_n | |
| 27 gamm_u, gamm_v | |
|
27
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
28 lambda |
| 0 | 29 end |
| 30 | |
| 31 methods | |
| 32 function obj = Wave2dCurve(m,ti,order,c,opSet) | |
| 33 default_arg('opSet',@sbp.Variable); | |
|
27
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
34 default_arg('c', 1); |
| 0 | 35 |
| 36 if length(m) == 1 | |
| 37 m = [m m]; | |
| 38 end | |
| 39 | |
| 40 m_u = m(1); | |
| 41 m_v = m(2); | |
| 42 m_tot = m_u*m_v; | |
| 43 | |
| 44 [u, h_u] = util.get_grid(0, 1, m_u); | |
| 45 [v, h_v] = util.get_grid(0, 1, m_v); | |
| 46 | |
| 47 | |
| 48 % Operators | |
| 49 ops_u = opSet(m_u,h_u,order); | |
| 50 ops_v = opSet(m_v,h_v,order); | |
| 51 | |
| 52 I_u = speye(m_u); | |
| 53 I_v = speye(m_v); | |
| 54 | |
| 55 D1_u = sparse(ops_u.derivatives.D1); | |
| 56 D2_u = ops_u.derivatives.D2; | |
| 57 H_u = sparse(ops_u.norms.H); | |
| 58 Hi_u = sparse(ops_u.norms.HI); | |
| 59 % M_u = sparse(ops_u.norms.M); | |
| 60 e_l_u = sparse(ops_u.boundary.e_1); | |
| 61 e_r_u = sparse(ops_u.boundary.e_m); | |
| 62 d1_l_u = sparse(ops_u.boundary.S_1); | |
| 63 d1_r_u = sparse(ops_u.boundary.S_m); | |
| 64 | |
| 65 D1_v = sparse(ops_v.derivatives.D1); | |
| 66 D2_v = ops_v.derivatives.D2; | |
| 67 H_v = sparse(ops_v.norms.H); | |
| 68 Hi_v = sparse(ops_v.norms.HI); | |
| 69 % M_v = sparse(ops_v.norms.M); | |
| 70 e_l_v = sparse(ops_v.boundary.e_1); | |
| 71 e_r_v = sparse(ops_v.boundary.e_m); | |
| 72 d1_l_v = sparse(ops_v.boundary.S_1); | |
| 73 d1_r_v = sparse(ops_v.boundary.S_m); | |
| 74 | |
| 75 | |
| 76 % Metric derivatives | |
| 77 [X,Y] = ti.map(u,v); | |
| 78 | |
| 79 [x_u,x_v] = gridDerivatives(X,D1_u,D1_v); | |
| 80 [y_u,y_v] = gridDerivatives(Y,D1_u,D1_v); | |
| 81 | |
| 82 | |
| 83 | |
| 84 J = x_u.*y_v - x_v.*y_u; | |
| 85 a11 = 1./J .* (x_v.^2 + y_v.^2); %% GÖR SOM MATRISER | |
| 86 a12 = -1./J .* (x_u.*x_v + y_u.*y_v); | |
| 87 a22 = 1./J .* (x_u.^2 + y_u.^2); | |
| 88 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2)); | |
| 89 | |
| 90 dof_order = reshape(1:m_u*m_v,m_v,m_u); | |
| 91 | |
| 92 Duu = sparse(m_tot); | |
| 93 Dvv = sparse(m_tot); | |
| 94 | |
| 95 for i = 1:m_v | |
| 96 D = D2_u(a11(i,:)); | |
| 97 p = dof_order(i,:); | |
| 98 Duu(p,p) = D; | |
| 99 end | |
| 100 | |
| 101 for i = 1:m_u | |
| 102 D = D2_v(a22(:,i)); | |
| 103 p = dof_order(:,i); | |
| 104 Dvv(p,p) = D; | |
| 105 end | |
| 106 | |
| 107 L_12 = spdiags(a12(:),0,m_tot,m_tot); | |
| 108 Du = kr(D1_u,I_v); | |
| 109 Dv = kr(I_u,D1_v); | |
| 110 | |
| 111 Duv = Du*L_12*Dv; | |
| 112 Dvu = Dv*L_12*Du; | |
| 113 | |
| 114 | |
| 115 | |
| 116 obj.H = kr(H_u,H_v); | |
| 117 obj.Hi = kr(Hi_u,Hi_v); | |
| 118 obj.Hu = kr(H_u,I_v); | |
| 119 obj.Hv = kr(I_u,H_v); | |
| 120 obj.Hiu = kr(Hi_u,I_v); | |
| 121 obj.Hiv = kr(I_u,Hi_v); | |
| 122 | |
| 123 % obj.M = kr(M_u,H_v)+kr(H_u,M_v); | |
| 124 obj.e_w = kr(e_l_u,I_v); | |
| 125 obj.e_e = kr(e_r_u,I_v); | |
| 126 obj.e_s = kr(I_u,e_l_v); | |
| 127 obj.e_n = kr(I_u,e_r_v); | |
| 128 obj.du_w = kr(d1_l_u,I_v); | |
| 129 obj.dv_w = (obj.e_w'*Dv)'; | |
| 130 obj.du_e = kr(d1_r_u,I_v); | |
| 131 obj.dv_e = (obj.e_e'*Dv)'; | |
| 132 obj.du_s = (obj.e_s'*Du)'; | |
| 133 obj.dv_s = kr(I_u,d1_l_v); | |
| 134 obj.du_n = (obj.e_n'*Du)'; | |
| 135 obj.dv_n = kr(I_u,d1_r_v); | |
| 136 | |
| 137 obj.m = m; | |
| 138 obj.h = [h_u h_v]; | |
| 139 obj.order = order; | |
| 140 | |
| 141 | |
| 142 obj.c = c; | |
| 143 obj.J = spdiags(J(:),0,m_tot,m_tot); | |
| 144 obj.Ji = spdiags(1./J(:),0,m_tot,m_tot); | |
| 145 obj.a11 = a11; | |
| 146 obj.a12 = a12; | |
| 147 obj.a22 = a22; | |
| 148 obj.D = obj.Ji*c^2*(Duu + Duv + Dvu + Dvv); | |
| 149 obj.u = u; | |
| 150 obj.v = v; | |
| 151 obj.X = X; | |
| 152 obj.Y = Y; | |
| 153 obj.x = X(:); | |
| 154 obj.y = Y(:); | |
|
27
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
155 obj.lambda = lambda; |
| 0 | 156 |
| 157 obj.gamm_u = h_u*ops_u.borrowing.M.S; | |
| 158 obj.gamm_v = h_v*ops_v.borrowing.M.S; | |
| 159 end | |
| 160 | |
| 161 | |
| 162 % Closure functions return the opertors applied to the own doamin to close the boundary | |
| 163 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. | |
| 164 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. | |
| 165 % type is a string specifying the type of boundary condition if there are several. | |
| 166 % data is a function returning the data that should be applied at the boundary. | |
| 167 % neighbour_scheme is an instance of Scheme that should be interfaced to. | |
| 168 % neighbour_boundary is a string specifying which boundary to interface to. | |
| 169 function [closure, penalty] = boundary_condition(obj,boundary,type,data) | |
| 170 default_arg('type','neumann'); | |
| 171 default_arg('data',0); | |
| 172 | |
| 173 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv] = obj.get_boundary_ops(boundary); | |
| 174 | |
| 175 switch type | |
| 176 % Dirichlet boundary condition | |
| 177 case {'D','d','dirichlet'} | |
| 178 error('not implemented') | |
| 179 alpha = obj.alpha; | |
| 180 | |
| 181 % tau1 < -alpha^2/gamma | |
| 182 tuning = 1.1; | |
| 183 tau1 = -tuning*alpha/gamm; | |
| 184 tau2 = s*alpha; | |
| 185 | |
| 186 p = tau1*e + tau2*d; | |
| 187 | |
| 188 closure = halfnorm_inv*p*e'; | |
| 189 | |
| 190 pp = halfnorm_inv*p; | |
| 191 switch class(data) | |
| 192 case 'double' | |
| 193 penalty = pp*data; | |
| 194 case 'function_handle' | |
| 195 penalty = @(t)pp*data(t); | |
| 196 otherwise | |
| 197 error('Weird data argument!') | |
| 198 end | |
| 199 | |
| 200 | |
| 201 % Neumann boundary condition | |
| 202 case {'N','n','neumann'} | |
| 203 c = obj.c; | |
| 204 | |
| 205 | |
| 206 a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); | |
| 207 a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); | |
| 208 d = (a_n * d_n' + a_t*d_t')'; | |
| 209 | |
| 210 tau1 = -s; | |
| 211 tau2 = 0; | |
| 212 tau = c.^2 * obj.Ji*(tau1*e + tau2*d); | |
| 213 | |
| 214 closure = halfnorm_inv*tau*d'; | |
| 215 | |
| 216 pp = halfnorm_inv*tau; | |
| 217 switch class(data) | |
| 218 case 'double' | |
| 219 penalty = pp*data; | |
| 220 case 'function_handle' | |
| 221 penalty = @(t)pp*data(t); | |
| 222 otherwise | |
| 223 error('Weird data argument!') | |
| 224 end | |
| 225 | |
| 226 % Unknown, boundary condition | |
| 227 otherwise | |
| 228 error('No such boundary condition: type = %s',type); | |
| 229 end | |
| 230 end | |
| 231 | |
| 232 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) | |
| 233 % u denotes the solution in the own domain | |
| 234 % v denotes the solution in the neighbour domain | |
|
5
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
235 tuning = 1.2; |
|
27
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
236 % tuning = 20.2; |
|
55
a8ed986fcf57
Minor renaming and clean up in 2d wave schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
27
diff
changeset
|
237 [e_u, d_n_u, d_t_u, coeff_n_u, coeff_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t] = obj.get_boundary_ops(boundary); |
|
a8ed986fcf57
Minor renaming and clean up in 2d wave schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
27
diff
changeset
|
238 [e_v, d_n_v, d_t_v, coeff_n_v, coeff_t_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t] = neighbour_scheme.get_boundary_ops(boundary); |
| 0 | 239 |
|
27
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
240 a_n_u = spdiag(coeff_n_u); |
|
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
241 a_t_u = spdiag(coeff_t_u); |
|
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
242 a_n_v = spdiag(coeff_n_v); |
|
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
243 a_t_v = spdiag(coeff_t_v); |
|
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
244 |
|
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
245 F_u = (s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u')'; |
|
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
246 F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')'; |
| 0 | 247 |
|
5
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
248 u = obj; |
|
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
249 v = neighbour_scheme; |
|
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
250 |
|
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
251 b1_u = gamm_u*u.lambda./u.a11.^2; |
|
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
252 b2_u = gamm_u*u.lambda./u.a22.^2; |
|
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
253 b1_v = gamm_v*v.lambda./v.a11.^2; |
|
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
254 b2_v = gamm_v*v.lambda./v.a22.^2; |
|
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
255 |
|
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
256 |
|
27
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
257 tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); |
|
5
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
258 m_tot = obj.m(1)*obj.m(2); |
|
27
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
259 tau = tuning * spdiag(tau(:)); |
| 0 | 260 sig1 = 1/2; |
|
27
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
261 sig2 = -1/2*s_u; |
| 0 | 262 |
|
27
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
263 penalty_parameter_1 = halfnorm_inv_u_n*(tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t)*e_u; |
| 0 | 264 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u; |
| 265 | |
| 266 | |
|
5
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
267 closure = obj.Ji*obj.c^2 * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u'); |
|
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
268 penalty = obj.Ji*obj.c^2 * (-penalty_parameter_1*e_v' - penalty_parameter_2*F_v'); |
| 0 | 269 end |
| 270 | |
| 271 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. | |
| 272 % The right boundary is considered the positive boundary | |
|
27
97a638f91fb8
Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
5
diff
changeset
|
273 function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,a_n, a_t] = get_boundary_ops(obj,boundary) |
| 0 | 274 switch boundary |
| 275 case 'w' | |
| 276 e = obj.e_w; | |
| 277 d_n = obj.du_w; | |
| 278 d_t = obj.dv_w; | |
| 279 s = -1; | |
| 280 | |
| 281 coeff_n = obj.a11(:,1); | |
| 282 coeff_t = obj.a12(:,1); | |
| 283 case 'e' | |
| 284 e = obj.e_e; | |
| 285 d_n = obj.du_e; | |
| 286 d_t = obj.dv_e; | |
| 287 s = 1; | |
| 288 | |
| 289 coeff_n = obj.a11(:,end); | |
| 290 coeff_t = obj.a12(:,end); | |
| 291 case 's' | |
| 292 e = obj.e_s; | |
| 293 d_n = obj.dv_s; | |
| 294 d_t = obj.du_s; | |
| 295 s = -1; | |
| 296 | |
| 297 coeff_n = obj.a22(1,:)'; | |
| 298 coeff_t = obj.a12(1,:)'; | |
| 299 case 'n' | |
| 300 e = obj.e_n; | |
| 301 d_n = obj.dv_n; | |
| 302 d_t = obj.du_n; | |
| 303 s = 1; | |
| 304 | |
| 305 coeff_n = obj.a22(end,:)'; | |
| 306 coeff_t = obj.a12(end,:)'; | |
| 307 otherwise | |
| 308 error('No such boundary: boundary = %s',boundary); | |
| 309 end | |
| 310 | |
| 311 switch boundary | |
| 312 case {'w','e'} | |
| 313 halfnorm_inv_n = obj.Hiu; | |
| 314 halfnorm_inv_t = obj.Hiv; | |
| 315 halfnorm_t = obj.Hv; | |
| 316 gamm = obj.gamm_u; | |
| 317 case {'s','n'} | |
| 318 halfnorm_inv_n = obj.Hiv; | |
| 319 halfnorm_inv_t = obj.Hiu; | |
| 320 halfnorm_t = obj.Hu; | |
| 321 gamm = obj.gamm_v; | |
| 322 end | |
| 323 end | |
| 324 | |
| 325 function N = size(obj) | |
| 326 N = prod(obj.m); | |
| 327 end | |
| 328 | |
| 329 end | |
| 330 | |
| 331 methods(Static) | |
| 332 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u | |
| 333 % and bound_v of scheme schm_v. | |
| 334 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') | |
| 335 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) | |
| 336 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); | |
| 337 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); | |
| 338 end | |
| 339 end | |
| 340 end |