Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2dCurve.m @ 298:861972361f75 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
The curvelinear works for quads
author Ylva Rydin <ylva.rydin@telia.com>
date Tue, 04 Oct 2016 08:40:42 +0200
parents
children 4d8d6eb0c116
comparison
equal deleted inserted replaced
297:cd30b22cee56 298:861972361f75
1 classdef Hypsyst2dCurve < scheme.Scheme
2 properties
3 m % Number of points in each direction, possibly a vector
4 n %size of system
5 h % Grid spacing
6 X,Y % Values of x and y for each grid point
7
8 J, Ji
9 xi,eta
10 Xi,Eta
11
12 A,B
13 X_eta, Y_eta
14 X_xi,Y_xi
15 order % Order accuracy for the approximation
16
17 D % non-stabalized scheme operator
18 Ahat, Bhat, E
19
20 H % Discrete norm
21 % Norms in the x and y directions
22 Hxii,Hetai % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
23 I_xi,I_eta, I_N, onesN
24 e_w, e_e, e_s, e_n
25 index_w, index_e,index_s,index_n
26 params %parameters for the coeficient matrice
27 end
28
29
30 methods
31 function obj = Hypsyst2dCurve(m, order, A, B, E, params, ti)
32 default_arg('E', [])
33 xilim = {0 1};
34 etalim = {0 1};
35
36 if length(m) == 1
37 m = [m m];
38 end
39 obj.params = params;
40 obj.A=A;
41 obj.B=B;
42
43
44 obj.Ahat=@(params,x,y,x_eta,y_eta)(A(params,x,y).*y_eta-B(params,x,y).*x_eta);
45 obj.Bhat=@(params,x,y,x_xi,y_xi)(B(params,x,y).*x_xi-A(params,x,y).*y_xi);
46 obj.E=@(params,x,y,~,~)E(params,x,y);
47
48 m_xi = m(1);
49 m_eta = m(2);
50 m_tot=m_xi*m_eta;
51
52 ops_xi = sbp.D2Standard(m_xi,xilim,order);
53 ops_eta = sbp.D2Standard(m_eta,etalim,order);
54
55 obj.xi = ops_xi.x;
56 obj.eta = ops_eta.x;
57
58 obj.Xi = kr(obj.xi,ones(m_eta,1));
59 obj.Eta = kr(ones(m_xi,1),obj.eta);
60
61 obj.n = length(A(obj.params,0,0));
62 obj.onesN=ones(obj.n);
63
64 obj.index_w=1:m_xi;
65 obj.index_e=(m_tot-m_xi+1):m_tot;
66 obj.index_s=1:m_xi:(m_tot-m_xi+1);
67 obj.index_n=(m_xi):m_xi:m_tot;
68
69 I_n = eye(obj.n);
70 I_xi = speye(m_xi);
71 obj.I_xi = I_xi;
72 I_eta = speye(m_eta);
73 obj.I_eta = I_eta;
74
75 D1_xi = kr(I_n, ops_xi.D1, I_eta);
76 obj.Hxii = kr(I_n, ops_xi.HI, I_eta);
77 D1_eta = kr(I_n, I_xi, ops_eta.D1);
78 obj.Hetai = kr(I_n, I_xi, ops_eta.HI);
79
80 obj.e_w = kr(I_n, ops_xi.e_l, I_eta);
81 obj.e_e = kr(I_n, ops_xi.e_r, I_eta);
82 obj.e_s = kr(I_n, I_xi, ops_eta.e_l);
83 obj.e_n = kr(I_n, I_xi, ops_eta.e_r);
84
85 [X,Y] = ti.map(obj.xi,obj.eta);
86
87 [x_xi,x_eta] = gridDerivatives(X,ops_xi.D1,ops_eta.D1);
88 [y_xi,y_eta] = gridDerivatives(Y,ops_xi.D1, ops_eta.D1);
89
90 obj.X=reshape(X,m_xi*m_eta,1);
91 obj.Y=reshape(Y,m_xi*m_eta,1);
92 obj.X_xi=reshape(x_xi,m_xi*m_eta,1);
93 obj.Y_xi=reshape(y_xi,m_xi*m_eta,1);
94 obj.X_eta=reshape(x_eta,m_xi*m_eta,1);
95 obj.Y_eta=reshape(y_eta,m_xi*m_eta,1);
96
97 Ahat_evaluated = obj.evaluateCoefficientMatrix(obj.Ahat, obj.X, obj.Y,obj.X_eta,obj.Y_eta);
98 Bhat_evaluated = obj.evaluateCoefficientMatrix(obj.Bhat, obj.X, obj.Y,obj.X_xi,obj.Y_xi);
99 E_evaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,[],[]);
100
101 obj.m=m;
102 obj.h=[ops_xi.h ops_eta.h];
103 obj.order=order;
104 obj.J=x_xi.*y_eta - x_eta.*y_xi;
105 obj.Ji =kr(I_n,spdiags(1./obj.J(:),0,m_xi*m_eta,m_xi*m_eta));
106
107 obj.D=obj.Ji*(-Ahat_evaluated*D1_xi-Bhat_evaluated*D1_eta)-E_evaluated;
108
109 end
110
111 % Closure functions return the opertors applied to the own doamin to close the boundary
112 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
113 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
114 % type is a string specifying the type of boundary condition if there are several.
115 % data is a function returning the data that should be applied at the boundary.
116 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
117 default_arg('type','char');
118 switch type
119 case{'c','char'}
120 [closure,penalty]=boundary_condition_char(obj,boundary);
121 case{'general'}
122 [closure,penalty]=boundary_condition_general(obj,boundary,L);
123 otherwise
124 error('No such boundary condition')
125 end
126 end
127
128 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
129 error('An interface function does not exist yet');
130 end
131
132 function N = size(obj)
133 N = obj.m;
134 end
135
136 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_)
137 params=obj.params;
138
139 if isa(mat,'function_handle')
140 [rows,cols]=size(mat(params,0,0,0,0));
141 x_=kr(x_,obj.onesN);
142 y_=kr(y_,obj.onesN);
143 matVec=mat(params,X',Y',x_',y_');
144 matVec=sparse(matVec);
145 side=max(length(X),length(Y));
146 else
147 matVec=mat;
148 [rows,cols]=size(matVec);
149 side=max(length(X),length(Y));
150 cols=cols/side;
151 end
152 ret=kron(ones(rows,cols),speye(side));
153
154 for ii=1:rows
155 for jj=1:cols
156 ret((ii-1)*side+1:ii*side,(jj-1)*side+1:jj*side)=diag(matVec(ii,(jj-1)*side+1:jj*side));
157 end
158 end
159 end
160
161
162 function [closure, penalty]=boundary_condition_char(obj,boundary)
163 params=obj.params;
164 xi=obj.xi; eta=obj.eta;
165 side=max(length(xi),length(eta));
166
167 switch boundary
168 case {'w','W','west'}
169 e_=obj.e_w;
170 mat=obj.Ahat;
171 boundPos='l';
172 Hi=obj.Hxii;
173 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi(1),eta,obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w));
174 case {'e','E','east'}
175 e_=obj.e_e;
176 mat=obj.Ahat;
177 boundPos='r';
178 Hi=obj.Hxii;
179 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi(end),eta,obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e));
180 case {'s','S','south'}
181 e_=obj.e_s;
182 mat=obj.Bhat;
183 boundPos='l';
184 Hi=obj.Hetai;
185 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi,eta(1),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s));
186 case {'n','N','north'}
187 e_=obj.e_n;
188 mat=obj.Bhat;
189 boundPos='r';
190 Hi=obj.Hetai;
191 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi,eta(end),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n));
192 end
193
194 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
195
196 switch boundPos
197 case {'l'}
198 tau=sparse(obj.n*side,pos*side);
199 Vi_plus=Vi(1:pos*side,:);
200 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side));
201 closure=Hi*e_*V*tau*Vi_plus*e_';
202 penalty=-Hi*e_*V*tau*Vi_plus;
203 case {'r'}
204 tau=sparse(obj.n*side,neg*side);
205 tau((pos+zeroval)*side+1:obj.n*side,:)=-abs(D((pos+zeroval)*side+1:obj.n*side,(pos+zeroval)*side+1:obj.n*side));
206 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:);
207 closure=Hi*e_*V*tau*Vi_minus*e_';
208 penalty=-Hi*e_*V*tau*Vi_minus;
209 end
210 end
211
212
213 function [closure,penalty]=boundary_condition_general(obj,boundary,L)
214 params=obj.params;
215 xi=obj.xi; eta=obj.eta;
216 side=max(length(xi),length(eta));
217
218 switch boundary
219 case {'w','W','west'}
220 e_=obj.e_w;
221 mat=obj.Ahat;
222 boundPos='l';
223 Hi=obj.Hxii;
224 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi(1),eta,obj.x_eta,obj.y_eta);
225 L=obj.evaluateCoefficientMatrix(L,xi(1),eta);
226 case {'e','E','east'}
227 e_=obj.e_e;
228 mat=obj.Ahat;
229 boundPos='r';
230 Hi=obj.Hxii;
231 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi(end),eta,obj.x_eta,obj.y_eta);
232 L=obj.evaluateCoefficientMatrix(L,xi(end),eta,[],[]);
233 case {'s','S','south'}
234 e_=obj.e_s;
235 mat=obj.Bhat;
236 boundPos='l';
237 Hi=obj.Hetai;
238 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi,eta(1),obj.x_xi,obj.y_xi);
239 L=obj.evaluateCoefficientMatrix(L,xi,eta(1));
240 case {'n','N','north'}
241 e_=obj.e_n;
242 mat=obj.Bhat;
243 boundPos='r';
244 Hi=obj.Hetai;
245 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi,eta(end),obj.x_xi,obj.y_xi);
246 L=obj.evaluateCoefficientMatrix(L,xi,eta(end));
247 end
248
249 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
250
251 switch boundPos
252 case {'l'}
253 tau=sparse(obj.n*side,pos*side);
254 Vi_plus=Vi(1:pos*side,:);
255 Vi_minus=Vi(pos*side+1:obj.n*side,:);
256 V_plus=V(:,1:pos*side);
257 V_minus=V(:,(pos+zeroval)*side+1:obj.n*side);
258
259 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side));
260 R=-inv(L*V_plus)*(L*V_minus);
261 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
262 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L;
263 case {'r'}
264 tau=sparse(obj.n*side,neg*side);
265 tau((pos+zeroval)*side+1:obj.n*side,:)=-abs(D((pos+zeroval)*side+1:obj.n*side,(pos+zeroval)*side+1:obj.n*side));
266 Vi_plus=Vi(1:pos*side,:);
267 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:);
268
269 V_plus=V(:,1:pos*side);
270 V_minus=V(:,(pos+zeroval)*side+1:obj.n*side);
271 R=-inv(L*V_minus)*(L*V_plus);
272 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
273 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L;
274 end
275 end
276
277
278 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,x_,y_)
279 params=obj.params;
280 syms xs ys
281 if(sum(abs(x_))~=0)
282 syms xs_
283 else
284 xs_=0;
285 end
286
287 if(sum(abs(y_))~=0)
288 syms ys_;
289 else
290 ys_=0;
291 end
292
293 [V, D]=eig(mat(params,xs,ys,xs_,ys_));
294 xs=1;ys=1; xs_=x_(1); ys_=y_(1);
295 DD=eval(diag(D));
296
297 poseig=find(DD>0);
298 zeroeig=find(DD==0);
299 negeig=find(DD<0);
300 syms xs ys xs_ ys_
301 DD=diag(D);
302
303 D=diag([DD(poseig);DD(zeroeig); DD(negeig)]);
304 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)];
305 Vi=inv(V);
306 xs=x;
307 ys=y;
308 xs_=x_';
309 ys_=y_';
310
311 side=max(length(x),length(y));
312 Dret=zeros(obj.n,side*obj.n);
313 Vret=zeros(obj.n,side*obj.n);
314 Viret=zeros(obj.n,side*obj.n);
315 for ii=1:obj.n
316 for jj=1:obj.n
317 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii));
318 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii));
319 Viret(jj,(ii-1)*side+1:side*ii)=eval(Vi(jj,ii));
320 end
321 end
322
323 D=sparse(Dret);
324 V=sparse(Vret);
325 Vi=sparse(Viret);
326 V=obj.evaluateCoefficientMatrix(V,x,y,x_,y_);
327 D=obj.evaluateCoefficientMatrix(D,x,y,x_,y_);
328 Vi=obj.evaluateCoefficientMatrix(Vi,x,y,x_,y_);
329 signVec=[length(poseig),length(zeroeig),length(negeig)];
330 end
331
332 end
333 end