Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst3d.m @ 349:cd6a29ab3746 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
A 3D is added and an attempt to imlement 3D transfinit interpolation has been initialized
author Ylva Rydin <ylva.rydin@telia.com>
date Thu, 13 Oct 2016 09:34:30 +0200
parents
children 5d5652fe826a
comparison
equal deleted inserted replaced
301:d9860ebc3148 349:cd6a29ab3746
1 classdef Hypsyst3d < 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, z % Grid
7 X, Y, Z% Values of x and y for each grid point
8 Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces
9 order % Order accuracy for the approximation
10
11 D % non-stabalized scheme operator
12 A, B, C, E
13
14 H % Discrete norm
15 % Norms in the x, y and z directions
16 Hxi,Hyi, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
17 I_x,I_y, I_z, I_N
18 e_w, e_e, e_s, e_n, e_b, e_t
19 params %parameters for the coeficient matrice
20 end
21
22
23 methods
24 function obj = Hypsyst3d(m, lim, order, A, B,C, E, params)
25 default_arg('E', [])
26 xlim = lim{1};
27 ylim = lim{2};
28 zlim = lim{3};
29
30 if length(m) == 1
31 m = [m m m];
32 end
33
34 obj.A=A;
35 obj.B=B;
36 obj.C=C;
37 obj.E=E;
38 m_x = m(1);
39 m_y = m(2);
40 m_z=m(3);
41 obj.params = params;
42
43 ops_x = sbp.D2Standard(m_x,xlim,order);
44 ops_y = sbp.D2Standard(m_y,ylim,order);
45 ops_z = sbp.D2Standard(m_z,zlim,order);
46
47 obj.x = ops_x.x;
48 obj.y = ops_y.x;
49 obj.z = ops_z.x;
50
51 obj.X = kr(obj.x,ones(m_y,1),ones(m_z,1));%% Que pasa?
52 obj.Y = kr(ones(m_x,1),obj.y,ones(m_z,1));
53 obj.Z = kr(ones(m_x,1),ones(m_y,1),obj.z);
54
55 obj.Yx=kr(obj.y,ones(m_z,1));
56 obj.Zx=kr(ones(m_y,1),obj.z);
57
58 obj.Xy=kr(obj.x,ones(m_z,1));
59 obj.Zy=kr(ones(m_x,1),obj.z);
60
61 obj.Xz=kr(obj.x,ones(m_y,1));
62 obj.Yz=kr(ones(m_z,1),obj.y);
63
64 Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z);
65 Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z);
66 Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z);
67 Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z);
68
69 obj.n = length(A(obj.params,0,0,0));
70
71 I_n = eye(obj.n);
72 I_x = speye(m_x);
73 obj.I_x = I_x;
74 I_y = speye(m_y);
75 obj.I_y = I_y;
76 I_z = speye(m_z);
77 obj.I_z = I_z;
78
79
80 D1_x = kr(I_n, ops_x.D1, I_y,I_z);
81 obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z);
82 D1_y = kr(I_n, I_x, ops_y.D1,I_z);
83 obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z);
84 D1_z = kr(I_n, I_x, I_y,ops_z.D1);
85 obj.Hzi = kr(I_n, I_x,I_y, ops_y.HI);
86
87 obj.e_w = kr(I_n, ops_x.e_l, I_y,I_z);
88 obj.e_e = kr(I_n, ops_x.e_r, I_y,I_z);
89 obj.e_s = kr(I_n, I_x, ops_y.e_l,I_z);
90 obj.e_n = kr(I_n, I_x, ops_y.e_r,I_z);
91 obj.e_b = kr(I_n, I_x, I_y, ops_z.e_l);
92 obj.e_t = kr(I_n, I_x, I_y, ops_z.e_r);
93
94 obj.m=m;
95 obj.h=[ops_x.h ops_y.h ops_x.h];
96 obj.order=order;
97
98 obj.D=-Aevaluated*D1_x-Bevaluated*D1_y-Cevaluated*D1_z-Eevaluated;
99 end
100
101 % Closure functions return the opertors applied to the own doamin to close the boundary
102 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
103 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
104 % type is a string specifying the type of boundary condition if there are several.
105 % data is a function returning the data that should be applied at the boundary.
106 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
107 default_arg('type','char');
108 BM=boundary_matrices(obj,boundary);
109
110 switch type
111 case{'c','char'}
112 [closure,penalty]=boundary_condition_char(obj,BM);
113 case{'general'}
114 [closure,penalty]=boundary_condition_general(obj,BM,boundary,L);
115 otherwise
116 error('No such boundary condition')
117 end
118 end
119
120 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
121 error('An interface function does not exist yet');
122 end
123
124 function N = size(obj)
125 N = obj.m;
126 end
127
128 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y, Z)
129 params=obj.params;
130 side=max(length(X),length(Y));
131 if isa(mat,'function_handle')
132 [rows,cols]=size(mat(params,0,0,0));
133 matVec=mat(params,X',Y',Z');
134 matVec=sparse(matVec);
135 else
136 matVec=mat;
137 [rows,cols]=size(matVec);
138 side=max(length(X),length(Y));
139 cols=cols/side;
140 end
141 ret=kron(ones(rows,cols),speye(side));
142
143 for ii=1:rows
144 for jj=1:cols
145 ret((ii-1)*side+1:ii*side,(jj-1)*side+1:jj*side)=diag(matVec(ii,(jj-1)*side+1:jj*side));
146 end
147 end
148 end
149
150
151 function [BM]=boundary_matrices(obj,boundary)
152 params=obj.params;
153
154 switch boundary
155 case {'w','W','west'}
156 BM.e_=obj.e_w;
157 mat=obj.A;
158 BM.boundpos='l';
159 BM.Hi=obj.Hxi;
160 [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.X(1),obj.Yx,obj.Zx);
161 BM.side=length(obj.Yx);
162 case {'e','E','east'}
163 BM.e_=obj.e_e;
164 mat=obj.A;
165 BM.boundpos='r';
166 BM.Hi=obj.Hxi;
167 [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.X(end),obj.Yx,obj.Zx);
168 BM.side=length(obj.Yx);
169 case {'s','S','south'}
170 BM.e_=obj.e_s;
171 mat=obj.B;
172 BM.boundpos='l';
173 BM.Hi=obj.Hyi;
174 [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.Xy,obj.Y(1),obj.Zy);
175 BM.side=length(obj.Xy);
176 case {'n','N','north'}
177 BM.e_=obj.e_n;
178 mat=obj.B;
179 BM.boundpos='r';
180 BM.Hi=obj.Hyi;
181 [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.Xy,obj.Y(end),obj.Zy);
182 BM.side=length(obj.Xy);
183 case{'b','B','Bottom'}
184 BM.e_=obj.e_b;
185 mat=obj.C;
186 BM.boundpos='l';
187 BM.Hi=obj.Hzi;
188 [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(1));
189 BM.side=length(obj.Xz);
190 case{'t','T','Top'}
191 BM.e_=obj.e_t;
192 mat=obj.C;
193 BM.boundpos='r';
194 BM.Hi=obj.Hzi;
195 [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(end));
196 BM.side=length(obj.Xz);
197 end
198
199 BM.pos=signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3);
200 end
201
202
203 function [closure, penalty]=boundary_condition_char(obj,BM)
204 side = BM.side;
205 pos = BM.pos;
206 neg = BM.neg;
207 zeroval=BM.zeroval;
208 V = BM.V;
209 Vi = BM.Vi;
210 Hi=BM.Hi;
211 D=BM.D;
212 e_=BM.e_;
213
214 switch BM.boundpos
215 case {'l'}
216 tau=sparse(obj.n*side,pos);
217 Vi_plus=Vi(1:pos,:);
218 tau(1:pos,:)=-abs(D(1:pos,1:pos));
219 closure=Hi*e_*V*tau*Vi_plus*e_';
220 penalty=-Hi*e_*V*tau*Vi_plus;
221 case {'r'}
222 tau=sparse(obj.n*side,neg);
223 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
224 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:);
225 closure=Hi*e_*V*tau*Vi_minus*e_';
226 penalty=-Hi*e_*V*tau*Vi_minus;
227 end
228 end
229
230
231 function [closure,penalty]=boundary_condition_general(obj,BM,boundary,L)
232 side = BM.side;
233 pos = BM.pos;
234 neg = BM.neg;
235 zeroval=BM.zeroval;
236 V = BM.V;
237 Vi = BM.Vi;
238 Hi=BM.Hi;
239 D=BM.D;
240 e_=BM.e_;
241 switch boundary
242 case {'w','W','west'}
243 L=obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx);
244 case {'e','E','east'}
245 L=obj.evaluateCoefficientMatrix(L,obj.x(end),obj.Yx,obj.Zx);
246 case {'s','S','south'}
247 L=obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(1),obj.Zy);
248 case {'n','N','north'}
249 L=obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(end),obj.Zy);
250 case {'b','B','bottom'}
251 L=obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(1));
252 case {'t','T','top'}
253 L=obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(end));
254 end
255
256 switch BM.boundpos
257 case {'l'}
258 tau=sparse(obj.n*side,pos);
259 Vi_plus=Vi(1:pos,:);
260 Vi_minus=Vi(pos+zeroval+1:obj.n*side,:);
261 V_plus=V(:,1:pos);
262 V_minus=V(:,(pos+zeroval)+1:obj.n*side);
263
264 tau(1:pos,:)=-abs(D(1:pos,1:pos));
265 R=-inv(L*V_plus)*(L*V_minus);
266 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
267 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L;
268 case {'r'}
269 tau=sparse(obj.n*side,neg);
270 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
271 Vi_plus=Vi(1:pos,:);
272 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:);
273
274 V_plus=V(:,1:pos);
275 V_minus=V(:,(pos+zeroval)+1:obj.n*side);
276 R=-inv(L*V_minus)*(L*V_plus);
277 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
278 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L;
279 end
280 end
281
282
283 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,z)
284 params=obj.params;
285 syms xs ys zs
286 [V, D]=eig(mat(params,xs,ys,zs));
287 xs=x;
288 ys=y;
289 zs=z;
290
291
292 side=max(length(x),length(y));
293 Dret=zeros(obj.n,side*obj.n);
294 Vret=zeros(obj.n,side*obj.n);
295 for ii=1:obj.n
296 for jj=1:obj.n
297 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii));
298 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii));
299 end
300 end
301
302 D=sparse(Dret);
303 V=sparse(Vret);
304 V=obj.evaluateCoefficientMatrix(V,x,y,z);
305 D=obj.evaluateCoefficientMatrix(D,x,y,z);
306 DD=diag(D);
307
308 poseig=(DD>0);
309 zeroeig=(DD==0);
310 negeig=(DD<0);
311
312 D=diag([DD(poseig); DD(zeroeig); DD(negeig)]);
313 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)];
314 Vi=inv(V);
315 signVec=[sum(poseig),sum(zeroeig),sum(negeig)];
316 end
317 end
318 end