Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst3d.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 9567a9dd220d
comparison
equal deleted inserted replaced
354:dbac99d2c318 355:69b078cf8072
24 24
25 methods 25 methods
26 function obj = Hypsyst3d(m, lim, order, A, B,C, E, params,operator) 26 function obj = Hypsyst3d(m, lim, order, A, B,C, E, params,operator)
27 default_arg('E', []) 27 default_arg('E', [])
28 default_arg('operatpr',[]) 28 default_arg('operatpr',[])
29 xlim = lim{1}; 29 xlim = lim{1};
30 ylim = lim{2}; 30 ylim = lim{2};
31 zlim = lim{3}; 31 zlim = lim{3};
32 32
33 if length(m) == 1 33 if length(m) == 1
34 m = [m m m]; 34 m = [m m m];
76 obj.Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z); 76 obj.Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z);
77 obj.Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z); 77 obj.Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z);
78 78
79 obj.n = length(A(obj.params,0,0,0)); 79 obj.n = length(A(obj.params,0,0,0));
80 80
81 I_n = eye(obj.n); 81 I_n = speye(obj.n);
82 I_x = speye(m_x); 82 I_x = speye(m_x);
83 obj.I_x = I_x; 83 obj.I_x = I_x;
84 I_y = speye(m_y); 84 I_y = speye(m_y);
85 obj.I_y = I_y; 85 obj.I_y = I_y;
86 I_z = speye(m_z); 86 I_z = speye(m_z);
88 I_N=kr(I_n,I_x,I_y,I_z); 88 I_N=kr(I_n,I_x,I_y,I_z);
89 89
90 obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z); 90 obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z);
91 obj.Hx = ops_x.H; 91 obj.Hx = ops_x.H;
92 obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z); 92 obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z);
93 obj.Hy = ops_y.H 93 obj.Hy = ops_y.H;
94 obj.Hzi = kr(I_n, I_x,I_y, ops_z.HI); 94 obj.Hzi = kr(I_n, I_x,I_y, ops_z.HI);
95 obj.Hz = ops_z.H; 95 obj.Hz = ops_z.H;
96 96
97 obj.e_w = kr(I_n, ops_x.e_l, I_y,I_z); 97 obj.e_w = kr(I_n, ops_x.e_l, I_y,I_z);
98 obj.e_e = kr(I_n, ops_x.e_r, I_y,I_z); 98 obj.e_e = kr(I_n, ops_x.e_r, I_y,I_z);
111 alphaB = max(eig(B(params,obj.x(end),obj.y(end),obj.z(end)))); 111 alphaB = max(eig(B(params,obj.x(end),obj.y(end),obj.z(end))));
112 alphaC = max(eig(C(params,obj.x(end),obj.y(end),obj.z(end)))); 112 alphaC = max(eig(C(params,obj.x(end),obj.y(end),obj.z(end))));
113 113
114 Ap = (obj.Aevaluated+alphaA*I_N)/2; 114 Ap = (obj.Aevaluated+alphaA*I_N)/2;
115 Am = (obj.Aevaluated-alphaA*I_N)/2; 115 Am = (obj.Aevaluated-alphaA*I_N)/2;
116 Dpx = kr(I_n, ops_x.Dp, I_y,I_z);
117 Dmx = kr(I_n, ops_x.Dm, I_y,I_z);
118 obj.D=-Am*Dpx;
119 temp=Ap*Dmx;
120 obj.D=obj.D-temp;
121 clear Ap Am Dpx Dmx
122
116 Bp = (obj.Bevaluated+alphaB*I_N)/2; 123 Bp = (obj.Bevaluated+alphaB*I_N)/2;
117 Bm = (obj.Bevaluated-alphaB*I_N)/2; 124 Bm = (obj.Bevaluated-alphaB*I_N)/2;
118 Cp = (obj.Cevaluated+alphaC*I_N)/2;
119 Cm = (obj.Cevaluated-alphaC*I_N)/2;
120
121 Dpx = kr(I_n, ops_x.Dp, I_y,I_z);
122 Dmx = kr(I_n, ops_x.Dm, I_y,I_z);
123 Dpy = kr(I_n, I_x, ops_y.Dp,I_z); 125 Dpy = kr(I_n, I_x, ops_y.Dp,I_z);
124 Dmy = kr(I_n, I_x, ops_y.Dm,I_z); 126 Dmy = kr(I_n, I_x, ops_y.Dm,I_z);
127 temp=Bm*Dpy;
128 obj.D=obj.D-temp;
129 temp=Bp*Dmy;
130 obj.D=obj.D-temp;
131 clear Bp Bm Dpy Dmy
132
133
134 Cp = (obj.Cevaluated+alphaC*I_N)/2;
135 Cm = (obj.Cevaluated-alphaC*I_N)/2;
125 Dpz = kr(I_n, I_x, I_y,ops_z.Dp); 136 Dpz = kr(I_n, I_x, I_y,ops_z.Dp);
126 Dmz = kr(I_n, I_x, I_y,ops_z.Dm); 137 Dmz = kr(I_n, I_x, I_y,ops_z.Dm);
127 138
128 obj.D=-Am*Dpx-Ap*Dmx-Bm*Dpy-Bp*Dmy-Cm*Dpz-Cp*Dmz-obj.Eevaluated; 139 temp=Cm*Dpz;
140 obj.D=obj.D-temp;
141 temp=-Cp*Dmz;
142 obj.D=obj.D-temp;
143 clear Cp Cm Dpz Dmz
144
145 obj.D=obj.D-obj.Eevaluated;
146
129 otherwise 147 otherwise
130 D1_x = kr(I_n, ops_x.D1, I_y,I_z); 148 D1_x = kr(I_n, ops_x.D1, I_y,I_z);
131 D1_y = kr(I_n, I_x, ops_y.D1,I_z); 149 D1_y = kr(I_n, I_x, ops_y.D1,I_z);
132 D1_z = kr(I_n, I_x, I_y,ops_z.D1); 150 D1_z = kr(I_n, I_x, I_y,ops_z.D1);
133 obj.D=-obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated; 151 obj.D=-obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated;