Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2d.m @ 297:cd30b22cee56 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
author Ylva Rydin <ylva.rydin@telia.com>
date Mon, 03 Oct 2016 08:33:47 +0200
parents a6ae1b104391
children d9860ebc3148
comparison
equal deleted inserted replaced
296:a6ae1b104391 297:cd30b22cee56
7 X,Y % Values of x and y for each grid point 7 X,Y % Values of x and y for each grid point
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 A, B, E 11 A, B, E
12 12
13 H % Discrete norm 13 H % Discrete norm
14 % Norms in the x and y directions 14 % Norms in the x and y directions
15 Hxi,Hyi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. 15 Hxi,Hyi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
16 I_x,I_y, I_N 16 I_x,I_y, I_N
17 e_w, e_e, e_s, e_n 17 e_w, e_e, e_s, e_n
18 params %parameters for the coeficient matrices 18 params %parameters for the coeficient matrice
19 matrices
20 end 19 end
21 20
22 21
23 methods 22 methods
24 function obj = Hypsyst2d(m, lim, order, A, B, E, params) 23 function obj = Hypsyst2d(m, lim, order, A, B, E, params)
27 ylim = lim{2}; 26 ylim = lim{2};
28 27
29 if length(m) == 1 28 if length(m) == 1
30 m = [m m]; 29 m = [m m];
31 end 30 end
31
32 obj.A=A;
33 obj.B=B;
34 obj.E=E;
32 35
33 m_x = m(1); 36 m_x = m(1);
34 m_y = m(2); 37 m_y = m(2);
35 obj.params = params; 38 obj.params = params;
36 39
37 obj.matrices = matrices;
38
39 ops_x = sbp.D2Standard(m_x,xlim,order); 40 ops_x = sbp.D2Standard(m_x,xlim,order);
40 ops_y = sbp.D2Standard(m_y,ylim,order); 41 ops_y = sbp.D2Standard(m_y,ylim,order);
41 42
42 obj.x = ops_x.x; 43 obj.x = ops_x.x;
43 obj.y = ops_y.x; 44 obj.y = ops_y.x;
44 45
45 obj.X = kr(obj.x,ones(m_y,1)); 46 obj.X = kr(obj.x,ones(m_y,1));
46 obj.Y = kr(ones(m_x,1),obj.y); 47 obj.Y = kr(ones(m_x,1),obj.y);
47 48
48 obj.A = obj.evaluateCoefficientMatrix(matrices.A, obj.X, obj.Y); 49 Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y);
49 obj.B = obj.evaluateCoefficientMatrix(matrices.B, obj.X, obj.Y); 50 Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y);
50 obj.E = obj.evaluateCoefficientMatrix(matrices.E, obj.X, obj.Y); 51 Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y);
51 52
52 obj.n = length(matrices.A(obj.params,0,0)); 53 obj.n = length(A(obj.params,0,0));
53 54
54 I_n = eye(obj.n);I_x = speye(m_x); 55 I_n = eye(obj.n);I_x = speye(m_x);
55 obj.I_x = I_x; 56 obj.I_x = I_x;
56 I_y = speye(m_y); 57 I_y = speye(m_y);
57 obj.I_y = I_y; 58 obj.I_y = I_y;
58 59
59 60
60 D1_x = kr(I_n, ops_x.D1, I_y); 61 D1_x = kr(I_n, ops_x.D1, I_y);
61 obj.Hxi = kr(I_n, ops_x.HI, I_y); 62 obj.Hxi = kr(I_n, ops_x.HI, I_y);
62 D1_y = kr(I_n, I_x, ops_y.D1)); 63 D1_y = kr(I_n, I_x, ops_y.D1);
63 obj.Hyi = kr(I_n, I_x, ops_y.HI)); 64 obj.Hyi = kr(I_n, I_x, ops_y.HI);
64 65
65 obj.e_w = kr(I_n, ops_x.e_l, I_y); 66 obj.e_w = kr(I_n, ops_x.e_l, I_y);
66 obj.e_e = kr(I_n, ops_x.e_r, I_y); 67 obj.e_e = kr(I_n, ops_x.e_r, I_y);
67 obj.e_s = kr(I_n, I_x, ops_y.e_l); 68 obj.e_s = kr(I_n, I_x, ops_y.e_l);
68 obj.e_n = kr(I_n, I_x, ops_y.e_r); 69 obj.e_n = kr(I_n, I_x, ops_y.e_r);
69 70
70 obj.m=m; 71 obj.m=m;
71 obj.h=[ops_x.h ops_y.h]; 72 obj.h=[ops_x.h ops_y.h];
72 obj.order=order; 73 obj.order=order;
73 74
74 obj.D=-obj.A*D1_x-obj.B*D1_y-obj.E; 75 obj.D=-Aevaluated*D1_x-Bevaluated*D1_y-Eevaluated;
75 76
76 end 77 end
77 78
78 % Closure functions return the opertors applied to the own doamin to close the boundary 79 % Closure functions return the opertors applied to the own doamin to close the boundary
79 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 80 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
130 side=max(length(x),length(y)); 131 side=max(length(x),length(y));
131 132
132 switch boundary 133 switch boundary
133 case {'w','W','west'} 134 case {'w','W','west'}
134 e_=obj.e_w; 135 e_=obj.e_w;
135 mat=obj.matrices.A; 136 mat=obj.A;
136 boundPos='l'; 137 boundPos='l';
137 Hi=obj.Hxi; 138 Hi=obj.Hxi;
138 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y); 139 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y);
139 case {'e','E','east'} 140 case {'e','E','east'}
140 e_=obj.e_e; 141 e_=obj.e_e;
141 mat=obj.matrices.A; 142 mat=obj.A;
142 boundPos='r'; 143 boundPos='r';
143 Hi=obj.Hxi; 144 Hi=obj.Hxi;
144 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y); 145 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y);
145 case {'s','S','south'} 146 case {'s','S','south'}
146 e_=obj.e_s; 147 e_=obj.e_s;
147 mat=obj.matrices.B; 148 mat=obj.B;
148 boundPos='l'; 149 boundPos='l';
149 Hi=obj.Hxi; 150 Hi=obj.Hyi;
150 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1)); 151 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1));
151 case {'n','N','north'} 152 case {'n','N','north'}
152 e_=obj.e_n; 153 e_=obj.e_n;
153 mat=obj.matrices.B; 154 mat=obj.B;
154 boundPos='r'; 155 boundPos='r';
155 Hi=obj.Hxi; 156 Hi=obj.Hyi;
156 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end)); 157 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end));
157 end 158 end
158 159
159 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 160 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
160 161
181 side=max(length(x),length(y)); 182 side=max(length(x),length(y));
182 183
183 switch boundary 184 switch boundary
184 case {'w','W','west'} 185 case {'w','W','west'}
185 e_=obj.e_w; 186 e_=obj.e_w;
186 mat=obj.matrices.A; 187 mat=obj.A;
187 boundPos='l'; 188 boundPos='l';
188 Hi=obj.Hxi; 189 Hi=obj.Hxi;
189 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y); 190 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y);
190 L=obj.evaluateCoefficientMatrix(L,x(1),y); 191 L=obj.evaluateCoefficientMatrix(L,x(1),y);
191 case {'e','E','east'} 192 case {'e','E','east'}
192 e_=obj.e_e; 193 e_=obj.e_e;
193 mat=obj.matrices.A; 194 mat=obj.A;
194 boundPos='r'; 195 boundPos='r';
195 Hi=obj.Hxi; 196 Hi=obj.Hxi;
196 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y); 197 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y);
197 L=obj.evaluateCoefficientMatrix(L,x(end),y); 198 L=obj.evaluateCoefficientMatrix(L,x(end),y);
198 case {'s','S','south'} 199 case {'s','S','south'}
199 e_=obj.e_s; 200 e_=obj.e_s;
200 mat=obj.matrices.B; 201 mat=obj.B;
201 boundPos='l'; 202 boundPos='l';
202 Hi=obj.Hxi; 203 Hi=obj.Hyi;
203 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1)); 204 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1));
204 L=obj.evaluateCoefficientMatrix(L,x,y(1)); 205 L=obj.evaluateCoefficientMatrix(L,x,y(1));
205 case {'n','N','north'} 206 case {'n','N','north'}
206 e_=obj.e_n; 207 e_=obj.e_n;
207 mat=obj.matrices.B; 208 mat=obj.B;
208 boundPos='r'; 209 boundPos='r';
209 Hi=obj.Hxi; 210 Hi=obj.Hyi;
210 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end)); 211 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end));
211 L=obj.evaluateCoefficientMatrix(L,x,y(end)); 212 L=obj.evaluateCoefficientMatrix(L,x,y(end));
212 end 213 end
213 214
214 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 215 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);