Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst3dCurve.m @ 355:69b078cf8072 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Upwind doed not work in the non curve-linear case.
author Ylva Rydin <ylva.rydin@telia.com>
date Thu, 08 Dec 2016 10:42:11 +0100
parents dbac99d2c318
children 7ada2db63268
comparison
equal deleted inserted replaced
354:dbac99d2c318 355:69b078cf8072
51 m_tot = m_xi*m_eta*m_zeta; 51 m_tot = m_xi*m_eta*m_zeta;
52 obj.params = params; 52 obj.params = params;
53 obj.n = length(A(obj,0,0,0)); 53 obj.n = length(A(obj,0,0,0));
54 54
55 obj.m = m; 55 obj.m = m;
56
57 obj.order = order; 56 obj.order = order;
58 obj.onesN = ones(obj.n); 57 obj.onesN = ones(obj.n);
59 58
60 switch operator 59 switch operator
61 case 'upwind' 60 case 'upwind'
183 alphaA = max(eig(obj.Ahat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_eta(end),obj.X_zeta(end),obj.Y_eta(end),obj.Y_zeta(end),obj.Z_eta(end),obj.Z_zeta(end)))); 182 alphaA = max(eig(obj.Ahat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_eta(end),obj.X_zeta(end),obj.Y_eta(end),obj.Y_zeta(end),obj.Z_eta(end),obj.Z_zeta(end))));
184 alphaB = max(eig(obj.Bhat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_zeta(end),obj.X_xi(end),obj.Y_zeta(end),obj.Y_xi(end),obj.Z_zeta(end),obj.Z_xi(end)))); 183 alphaB = max(eig(obj.Bhat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_zeta(end),obj.X_xi(end),obj.Y_zeta(end),obj.Y_xi(end),obj.Z_zeta(end),obj.Z_xi(end))));
185 alphaC = max(eig(obj.Chat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_xi(end),obj.X_eta(end),obj.Y_xi(end),obj.Y_eta(end),obj.Z_xi(end),obj.Z_eta(end)))); 184 alphaC = max(eig(obj.Chat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_xi(end),obj.X_eta(end),obj.Y_xi(end),obj.Y_eta(end),obj.Z_xi(end),obj.Z_eta(end))));
186 185
187 Ap = (obj.Aevaluated+alphaA*I_N)/2; 186 Ap = (obj.Aevaluated+alphaA*I_N)/2;
187 Dmxi = kr(I_n, ops_xi.Dm, I_eta,I_zeta);
188 diffSum=-Ap*Dmxi;
189 clear Ap Dmxi
190
188 Am = (obj.Aevaluated-alphaA*I_N)/2; 191 Am = (obj.Aevaluated-alphaA*I_N)/2;
192 clear obj.Aevaluated
193 Dpxi = kr(I_n, ops_xi.Dp, I_eta,I_zeta);
194 temp=Am*Dpxi;
195 diffSum=diffSum-temp;
196 clear Am Dpxi
197
189 Bp = (obj.Bevaluated+alphaB*I_N)/2; 198 Bp = (obj.Bevaluated+alphaB*I_N)/2;
199 Dmeta = kr(I_n, I_xi, ops_eta.Dm,I_zeta);
200 temp=Bp*Dmeta;
201 diffSum=diffSum-temp;
202 clear Bp Dmeta
203
190 Bm = (obj.Bevaluated-alphaB*I_N)/2; 204 Bm = (obj.Bevaluated-alphaB*I_N)/2;
205 clear obj.Bevaluated
206 Dpeta = kr(I_n, I_xi, ops_eta.Dp,I_zeta);
207 temp=Bm*Dpeta;
208 diffSum=diffSum-temp;
209
210
191 Cp = (obj.Cevaluated+alphaC*I_N)/2; 211 Cp = (obj.Cevaluated+alphaC*I_N)/2;
192 Cm = (obj.Cevaluated-alphaC*I_N)/2; 212 Dmzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dm);
193 213 temp=Cp*Dmzeta;
194 Dpxi = kr(I_n, ops_xi.Dp, I_eta,I_zeta); 214 diffSum=diffSum-temp;
195 Dmxi = kr(I_n, ops_xi.Dm, I_eta,I_zeta); 215 clear Cp Dmzeta
196 Dpeta = kr(I_n, I_xi, ops_eta.Dp,I_zeta); 216
197 Dmeta = kr(I_n, I_xi, ops_eta.Dm,I_zeta); 217 Cm = (obj.Cevaluated-alphaC*I_N)/2;
218 clear obj.Cevaluated
198 Dpzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dp); 219 Dpzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dp);
199 Dmzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dm); 220 temp=Cm*Dpzeta;
200 221 diffSum=diffSum-temp;
201 obj.D = obj.Ji*(-Am*Dpxi-Ap*Dmxi-Bm*Dpeta-Bp*Dmeta-Cm*Dpzeta-Cp*Dmzeta)-obj.Eevaluated; 222
223 obj.D = obj.Ji*diffSum-obj.Eevaluated;
202 otherwise 224 otherwise
203 obj.D = obj.Ji*(-obj.Aevaluated*D1_xi-obj.Bevaluated*D1_eta -obj.Cevaluated*D1_zeta)-obj.Eevaluated; 225 obj.D = obj.Ji*(-obj.Aevaluated*D1_xi-obj.Bevaluated*D1_eta -obj.Cevaluated*D1_zeta)-obj.Eevaluated;
204 end 226 end
205 227
206 228