Mercurial > repos > public > sbplib
annotate +scheme/Wave2d.m @ 577:e45c9b56d50d feature/grids
Add an Empty grid class
The need turned up for the flexural code when we may or may not have a grid for the open water and want to plot that solution.
In case there is no open water we need an empty grid to plot the empty gridfunction against to avoid errors.
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Thu, 07 Sep 2017 09:16:12 +0200 |
parents | a8ed986fcf57 |
children | 459eeb99130f |
rev | line source |
---|---|
55
a8ed986fcf57
Minor renaming and clean up in 2d wave schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
1 classdef Wave2d < scheme.Scheme |
0 | 2 properties |
3 m % Number of points in each direction, possibly a vector | |
4 h % Grid spacing | |
5 x,y % Grid | |
6 X,Y % Values of x and y for each grid point | |
7 order % Order accuracy for the approximation | |
8 | |
9 D % non-stabalized scheme operator | |
10 M % Derivative norm | |
11 alpha | |
12 | |
13 H % Discrete norm | |
14 Hi | |
15 H_x, H_y % Norms in the x and y directions | |
16 Hx,Hy % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. | |
17 Hi_x, Hi_y | |
18 Hix, Hiy | |
19 e_w, e_e, e_s, e_n | |
20 d1_w, d1_e, d1_s, d1_n | |
21 gamm_x, gamm_y | |
22 end | |
23 | |
24 methods | |
55
a8ed986fcf57
Minor renaming and clean up in 2d wave schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
25 function obj = Wave2d(m,lim,order,alpha) |
a8ed986fcf57
Minor renaming and clean up in 2d wave schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
26 default_arg('alpha',1); |
a8ed986fcf57
Minor renaming and clean up in 2d wave schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
27 |
a8ed986fcf57
Minor renaming and clean up in 2d wave schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
28 xlim = lim{1}; |
a8ed986fcf57
Minor renaming and clean up in 2d wave schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
29 ylim = lim{2}; |
0 | 30 |
31 if length(m) == 1 | |
32 m = [m m]; | |
33 end | |
34 | |
35 m_x = m(1); | |
36 m_y = m(2); | |
37 | |
38 [x, h_x] = util.get_grid(xlim{:},m_x); | |
39 [y, h_y] = util.get_grid(ylim{:},m_y); | |
40 | |
41 ops_x = sbp.Ordinary(m_x,h_x,order); | |
42 ops_y = sbp.Ordinary(m_y,h_y,order); | |
43 | |
44 I_x = speye(m_x); | |
45 I_y = speye(m_y); | |
46 | |
47 D2_x = sparse(ops_x.derivatives.D2); | |
48 H_x = sparse(ops_x.norms.H); | |
49 Hi_x = sparse(ops_x.norms.HI); | |
50 M_x = sparse(ops_x.norms.M); | |
51 e_l_x = sparse(ops_x.boundary.e_1); | |
52 e_r_x = sparse(ops_x.boundary.e_m); | |
53 d1_l_x = sparse(ops_x.boundary.S_1); | |
54 d1_r_x = sparse(ops_x.boundary.S_m); | |
55 | |
56 D2_y = sparse(ops_y.derivatives.D2); | |
57 H_y = sparse(ops_y.norms.H); | |
58 Hi_y = sparse(ops_y.norms.HI); | |
59 M_y = sparse(ops_y.norms.M); | |
60 e_l_y = sparse(ops_y.boundary.e_1); | |
61 e_r_y = sparse(ops_y.boundary.e_m); | |
62 d1_l_y = sparse(ops_y.boundary.S_1); | |
63 d1_r_y = sparse(ops_y.boundary.S_m); | |
64 | |
65 D2 = kr(D2_x, I_y) + kr(I_x, D2_y); | |
66 obj.H = kr(H_x,H_y); | |
67 obj.Hx = kr(H_x,I_y); | |
68 obj.Hy = kr(I_x,H_y); | |
69 obj.Hix = kr(Hi_x,I_y); | |
70 obj.Hiy = kr(I_x,Hi_y); | |
71 obj.Hi = kr(Hi_x,Hi_y); | |
72 obj.M = kr(M_x,H_y)+kr(H_x,M_y); | |
73 obj.e_w = kr(e_l_x,I_y); | |
74 obj.e_e = kr(e_r_x,I_y); | |
75 obj.e_s = kr(I_x,e_l_y); | |
76 obj.e_n = kr(I_x,e_r_y); | |
77 obj.d1_w = kr(d1_l_x,I_y); | |
78 obj.d1_e = kr(d1_r_x,I_y); | |
79 obj.d1_s = kr(I_x,d1_l_y); | |
80 obj.d1_n = kr(I_x,d1_r_y); | |
81 | |
82 obj.m = m; | |
83 obj.h = [h_x h_y]; | |
84 obj.order = order; | |
85 | |
86 obj.alpha = alpha; | |
87 obj.D = alpha*D2; | |
88 obj.x = x; | |
89 obj.y = y; | |
90 obj.X = kr(x,ones(m_y,1)); | |
91 obj.Y = kr(ones(m_x,1),y); | |
92 | |
93 obj.gamm_x = h_x*ops_x.borrowing.M.S; | |
94 obj.gamm_y = h_y*ops_y.borrowing.M.S; | |
95 end | |
96 | |
97 | |
98 % Closure functions return the opertors applied to the own doamin to close the boundary | |
99 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. | |
100 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. | |
101 % type is a string specifying the type of boundary condition if there are several. | |
102 % data is a function returning the data that should be applied at the boundary. | |
103 % neighbour_scheme is an instance of Scheme that should be interfaced to. | |
104 % neighbour_boundary is a string specifying which boundary to interface to. | |
105 function [closure, penalty] = boundary_condition(obj,boundary,type,data) | |
106 default_arg('type','neumann'); | |
107 default_arg('data',0); | |
108 | |
109 [e,d,s,gamm,halfnorm_inv] = obj.get_boundary_ops(boundary); | |
110 | |
111 switch type | |
112 % Dirichlet boundary condition | |
113 case {'D','d','dirichlet'} | |
114 alpha = obj.alpha; | |
115 | |
116 % tau1 < -alpha^2/gamma | |
117 tuning = 1.1; | |
118 tau1 = -tuning*alpha/gamm; | |
119 tau2 = s*alpha; | |
120 | |
121 p = tau1*e + tau2*d; | |
122 | |
123 closure = halfnorm_inv*p*e'; | |
124 | |
125 pp = halfnorm_inv*p; | |
126 switch class(data) | |
127 case 'double' | |
128 penalty = pp*data; | |
129 case 'function_handle' | |
130 penalty = @(t)pp*data(t); | |
131 otherwise | |
132 error('Wierd data argument!') | |
133 end | |
134 | |
135 | |
136 % Neumann boundary condition | |
137 case {'N','n','neumann'} | |
138 alpha = obj.alpha; | |
139 tau1 = -s*alpha; | |
140 tau2 = 0; | |
141 tau = tau1*e + tau2*d; | |
142 | |
143 closure = halfnorm_inv*tau*d'; | |
144 | |
145 pp = halfnorm_inv*tau; | |
146 switch class(data) | |
147 case 'double' | |
148 penalty = pp*data; | |
149 case 'function_handle' | |
150 penalty = @(t)pp*data(t); | |
151 otherwise | |
152 error('Wierd data argument!') | |
153 end | |
154 | |
155 % Unknown, boundary condition | |
156 otherwise | |
157 error('No such boundary condition: type = %s',type); | |
158 end | |
159 end | |
160 | |
161 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) | |
162 % u denotes the solution in the own domain | |
163 % v denotes the solution in the neighbour domain | |
164 [e_u,d_u,s_u,gamm_u, halfnorm_inv] = obj.get_boundary_ops(boundary); | |
165 [e_v,d_v,s_v,gamm_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); | |
166 | |
167 tuning = 1.1; | |
168 | |
169 alpha_u = obj.alpha; | |
170 alpha_v = neighbour_scheme.alpha; | |
171 | |
172 % tau1 < -(alpha_u/gamm_u + alpha_v/gamm_v) | |
173 | |
174 tau1 = -(alpha_u/gamm_u + alpha_v/gamm_v) * tuning; | |
175 tau2 = s_u*1/2*alpha_u; | |
176 sig1 = s_u*(-1/2); | |
177 sig2 = 0; | |
178 | |
179 tau = tau1*e_u + tau2*d_u; | |
180 sig = sig1*e_u + sig2*d_u; | |
181 | |
182 closure = halfnorm_inv*( tau*e_u' + sig*alpha_u*d_u'); | |
183 penalty = halfnorm_inv*(-tau*e_v' - sig*alpha_v*d_v'); | |
184 end | |
185 | |
186 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. | |
187 % The right boundary is considered the positive boundary | |
188 function [e,d,s,gamm, halfnorm_inv] = get_boundary_ops(obj,boundary) | |
189 switch boundary | |
190 case 'w' | |
191 e = obj.e_w; | |
192 d = obj.d1_w; | |
193 s = -1; | |
194 gamm = obj.gamm_x; | |
195 halfnorm_inv = obj.Hix; | |
196 case 'e' | |
197 e = obj.e_e; | |
198 d = obj.d1_e; | |
199 s = 1; | |
200 gamm = obj.gamm_x; | |
201 halfnorm_inv = obj.Hix; | |
202 case 's' | |
203 e = obj.e_s; | |
204 d = obj.d1_s; | |
205 s = -1; | |
206 gamm = obj.gamm_y; | |
207 halfnorm_inv = obj.Hiy; | |
208 case 'n' | |
209 e = obj.e_n; | |
210 d = obj.d1_n; | |
211 s = 1; | |
212 gamm = obj.gamm_y; | |
213 halfnorm_inv = obj.Hiy; | |
214 otherwise | |
215 error('No such boundary: boundary = %s',boundary); | |
216 end | |
217 end | |
218 | |
219 function N = size(obj) | |
220 N = prod(obj.m); | |
221 end | |
222 | |
223 end | |
224 | |
225 methods(Static) | |
226 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u | |
227 % and bound_v of scheme schm_v. | |
228 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') | |
229 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) | |
230 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); | |
231 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); | |
232 end | |
233 end | |
234 end |