Mercurial > repos > public > sbplib
changeset 355:69b078cf8072 feature/hypsyst
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 | 9567a9dd220d |
| files | +grid/Ti3D.m +scheme/Hypsyst3d.m +scheme/Hypsyst3dCurve.m |
| diffstat | 3 files changed, 63 insertions(+), 22 deletions(-) [+] |
line wrap: on
line diff
diff -r dbac99d2c318 -r 69b078cf8072 +grid/Ti3D.m --- a/+grid/Ti3D.m Mon Nov 28 08:46:28 2016 +0100 +++ b/+grid/Ti3D.m Thu Dec 08 10:42:11 2016 +0100 @@ -5,8 +5,9 @@ end methods - % TODO function to label boundary names. - % function to find largest and smallest delta h in the grid. Maybe shouldnt live here + % TODO write all fancy features for flipping around with the surfaces + % Each surface is defined with an outward facing outward and choosing + % the "corner" where XI=0 if not possible the corner where ETA=0 is choosen function obj = Ti3D(CW,CE,CS,CN,CB,CT) obj.gs = {CE,CW,CS,CN,CB,CT};
diff -r dbac99d2c318 -r 69b078cf8072 +scheme/Hypsyst3d.m --- a/+scheme/Hypsyst3d.m Mon Nov 28 08:46:28 2016 +0100 +++ b/+scheme/Hypsyst3d.m Thu Dec 08 10:42:11 2016 +0100 @@ -26,7 +26,7 @@ function obj = Hypsyst3d(m, lim, order, A, B,C, E, params,operator) default_arg('E', []) default_arg('operatpr',[]) - xlim = lim{1}; + xlim = lim{1}; ylim = lim{2}; zlim = lim{3}; @@ -78,7 +78,7 @@ obj.n = length(A(obj.params,0,0,0)); - I_n = eye(obj.n); + I_n = speye(obj.n); I_x = speye(m_x); obj.I_x = I_x; I_y = speye(m_y); @@ -90,7 +90,7 @@ obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z); obj.Hx = ops_x.H; obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z); - obj.Hy = ops_y.H + obj.Hy = ops_y.H; obj.Hzi = kr(I_n, I_x,I_y, ops_z.HI); obj.Hz = ops_z.H; @@ -113,19 +113,37 @@ Ap = (obj.Aevaluated+alphaA*I_N)/2; Am = (obj.Aevaluated-alphaA*I_N)/2; + Dpx = kr(I_n, ops_x.Dp, I_y,I_z); + Dmx = kr(I_n, ops_x.Dm, I_y,I_z); + obj.D=-Am*Dpx; + temp=Ap*Dmx; + obj.D=obj.D-temp; + clear Ap Am Dpx Dmx + Bp = (obj.Bevaluated+alphaB*I_N)/2; Bm = (obj.Bevaluated-alphaB*I_N)/2; - Cp = (obj.Cevaluated+alphaC*I_N)/2; - Cm = (obj.Cevaluated-alphaC*I_N)/2; - - Dpx = kr(I_n, ops_x.Dp, I_y,I_z); - Dmx = kr(I_n, ops_x.Dm, I_y,I_z); Dpy = kr(I_n, I_x, ops_y.Dp,I_z); Dmy = kr(I_n, I_x, ops_y.Dm,I_z); + temp=Bm*Dpy; + obj.D=obj.D-temp; + temp=Bp*Dmy; + obj.D=obj.D-temp; + clear Bp Bm Dpy Dmy + + + Cp = (obj.Cevaluated+alphaC*I_N)/2; + Cm = (obj.Cevaluated-alphaC*I_N)/2; Dpz = kr(I_n, I_x, I_y,ops_z.Dp); Dmz = kr(I_n, I_x, I_y,ops_z.Dm); - obj.D=-Am*Dpx-Ap*Dmx-Bm*Dpy-Bp*Dmy-Cm*Dpz-Cp*Dmz-obj.Eevaluated; + temp=Cm*Dpz; + obj.D=obj.D-temp; + temp=-Cp*Dmz; + obj.D=obj.D-temp; + clear Cp Cm Dpz Dmz + + obj.D=obj.D-obj.Eevaluated; + otherwise D1_x = kr(I_n, ops_x.D1, I_y,I_z); D1_y = kr(I_n, I_x, ops_y.D1,I_z);
diff -r dbac99d2c318 -r 69b078cf8072 +scheme/Hypsyst3dCurve.m --- a/+scheme/Hypsyst3dCurve.m Mon Nov 28 08:46:28 2016 +0100 +++ b/+scheme/Hypsyst3dCurve.m Thu Dec 08 10:42:11 2016 +0100 @@ -53,7 +53,6 @@ obj.n = length(A(obj,0,0,0)); obj.m = m; - obj.order = order; obj.onesN = ones(obj.n); @@ -185,20 +184,43 @@ 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)))); Ap = (obj.Aevaluated+alphaA*I_N)/2; + Dmxi = kr(I_n, ops_xi.Dm, I_eta,I_zeta); + diffSum=-Ap*Dmxi; + clear Ap Dmxi + Am = (obj.Aevaluated-alphaA*I_N)/2; - Bp = (obj.Bevaluated+alphaB*I_N)/2; - Bm = (obj.Bevaluated-alphaB*I_N)/2; - Cp = (obj.Cevaluated+alphaC*I_N)/2; - Cm = (obj.Cevaluated-alphaC*I_N)/2; - + clear obj.Aevaluated Dpxi = kr(I_n, ops_xi.Dp, I_eta,I_zeta); - Dmxi = kr(I_n, ops_xi.Dm, I_eta,I_zeta); + temp=Am*Dpxi; + diffSum=diffSum-temp; + clear Am Dpxi + + Bp = (obj.Bevaluated+alphaB*I_N)/2; + Dmeta = kr(I_n, I_xi, ops_eta.Dm,I_zeta); + temp=Bp*Dmeta; + diffSum=diffSum-temp; + clear Bp Dmeta + + Bm = (obj.Bevaluated-alphaB*I_N)/2; + clear obj.Bevaluated Dpeta = kr(I_n, I_xi, ops_eta.Dp,I_zeta); - Dmeta = kr(I_n, I_xi, ops_eta.Dm,I_zeta); + temp=Bm*Dpeta; + diffSum=diffSum-temp; + + + Cp = (obj.Cevaluated+alphaC*I_N)/2; + Dmzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dm); + temp=Cp*Dmzeta; + diffSum=diffSum-temp; + clear Cp Dmzeta + + Cm = (obj.Cevaluated-alphaC*I_N)/2; + clear obj.Cevaluated Dpzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dp); - Dmzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dm); - - obj.D = obj.Ji*(-Am*Dpxi-Ap*Dmxi-Bm*Dpeta-Bp*Dmeta-Cm*Dpzeta-Cp*Dmzeta)-obj.Eevaluated; + temp=Cm*Dpzeta; + diffSum=diffSum-temp; + + obj.D = obj.Ji*diffSum-obj.Eevaluated; otherwise obj.D = obj.Ji*(-obj.Aevaluated*D1_xi-obj.Bevaluated*D1_eta -obj.Cevaluated*D1_zeta)-obj.Eevaluated; end