Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 553:63d5c07943dd feature/grids/laplace_refactor

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Reorder and restructure constructor
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 29 Aug 2017 11:04:15 +0200
parents ffaf13533c27
children 6e71d4f8b155
comparison
equal deleted inserted replaced
552:ffaf13533c27 553:63d5c07943dd
47 function obj = LaplaceCurvilinear(g ,order, a, b, opSet) 47 function obj = LaplaceCurvilinear(g ,order, a, b, opSet)
48 default_arg('opSet',@sbp.D2Variable); 48 default_arg('opSet',@sbp.D2Variable);
49 default_arg('a', 1); 49 default_arg('a', 1);
50 default_arg('b', 1); 50 default_arg('b', 1);
51 51
52
53 if b ~=1 52 if b ~=1
54 error('Not implemented yet') 53 error('Not implemented yet')
55 end 54 end
56 55
57 assert(isa(g, 'grid.Curvilinear')) 56 assert(isa(g, 'grid.Curvilinear'))
63 62
64 h = g.scaling(); 63 h = g.scaling();
65 h_u = h(1); 64 h_u = h(1);
66 h_v = h(2); 65 h_v = h(2);
67 66
68 % Operators 67
68 % 1D operators
69 ops_u = opSet(m_u, {0, 1}, order); 69 ops_u = opSet(m_u, {0, 1}, order);
70 ops_v = opSet(m_v, {0, 1}, order); 70 ops_v = opSet(m_v, {0, 1}, order);
71 71
72 I_u = speye(m_u); 72 I_u = speye(m_u);
73 I_v = speye(m_v); 73 I_v = speye(m_v);
88 e_l_v = ops_v.e_l; 88 e_l_v = ops_v.e_l;
89 e_r_v = ops_v.e_r; 89 e_r_v = ops_v.e_r;
90 d1_l_v = ops_v.d1_l; 90 d1_l_v = ops_v.d1_l;
91 d1_r_v = ops_v.d1_r; 91 d1_r_v = ops_v.d1_r;
92 92
93
94 % Logical operators
93 Du = kr(D1_u,I_v); 95 Du = kr(D1_u,I_v);
94 Dv = kr(I_u,D1_v); 96 Dv = kr(I_u,D1_v);
95 97 obj.Hu = kr(H_u,I_v);
96 % Metric derivatives 98 obj.Hv = kr(I_u,H_v);
99 obj.Hiu = kr(Hi_u,I_v);
100 obj.Hiv = kr(I_u,Hi_v);
101
102 e_w = kr(e_l_u,I_v);
103 e_e = kr(e_r_u,I_v);
104 e_s = kr(I_u,e_l_v);
105 e_n = kr(I_u,e_r_v);
106 obj.du_w = kr(d1_l_u,I_v);
107 obj.dv_w = (e_w'*Dv)';
108 obj.du_e = kr(d1_r_u,I_v);
109 obj.dv_e = (e_e'*Dv)';
110 obj.du_s = (e_s'*Du)';
111 obj.dv_s = kr(I_u,d1_l_v);
112 obj.du_n = (e_n'*Du)';
113 obj.dv_n = kr(I_u,d1_r_v);
114
115
116 % Metric coefficients
97 coords = g.points(); 117 coords = g.points();
98 x = coords(:,1); 118 x = coords(:,1);
99 y = coords(:,2); 119 y = coords(:,2);
100 120
101 x_u = Du*x; 121 x_u = Du*x;
107 a11 = 1./J .* (x_v.^2 + y_v.^2); 127 a11 = 1./J .* (x_v.^2 + y_v.^2);
108 a12 = -1./J .* (x_u.*x_v + y_u.*y_v); 128 a12 = -1./J .* (x_u.*x_v + y_u.*y_v);
109 a22 = 1./J .* (x_u.^2 + y_u.^2); 129 a22 = 1./J .* (x_u.^2 + y_u.^2);
110 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2)); 130 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2));
111 131
132 obj.x_u = x_u;
133 obj.x_v = x_v;
134 obj.y_u = y_u;
135 obj.y_v = y_v;
136
137
112 % Assemble full operators 138 % Assemble full operators
113 L_12 = spdiag(a12); 139 L_12 = spdiag(a12);
114 Duv = Du*L_12*Dv; 140 Duv = Du*L_12*Dv;
115 Dvu = Dv*L_12*Du; 141 Dvu = Dv*L_12*Du;
116 142
128 D = D2_v(a22(ind(i,:))); 154 D = D2_v(a22(ind(i,:)));
129 p = ind(i,:); 155 p = ind(i,:);
130 Dvv(p,p) = D; 156 Dvv(p,p) = D;
131 end 157 end
132 158
133 obj.H = kr(H_u,H_v); 159
134 obj.Hi = kr(Hi_u,Hi_v); 160 % Physical operators
135 obj.Hu = kr(H_u,I_v); 161 obj.J = spdiag(J);
136 obj.Hv = kr(I_u,H_v); 162 obj.Ji = spdiag(1./J);
137 obj.Hiu = kr(Hi_u,I_v); 163
138 obj.Hiv = kr(I_u,Hi_v); 164 obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv);
139 165 obj.H = obj.J*kr(H_u,H_v);
140 obj.e_w = kr(e_l_u,I_v); 166 obj.Hi = obj.Ji*kr(Hi_u,Hi_v);
141 obj.e_e = kr(e_r_u,I_v); 167
142 obj.e_s = kr(I_u,e_l_v); 168 obj.e_w = e_w;
143 obj.e_n = kr(I_u,e_r_v); 169 obj.e_e = e_e;
144 obj.du_w = kr(d1_l_u,I_v); 170 obj.e_s = e_s;
145 obj.dv_w = (obj.e_w'*Dv)'; 171 obj.e_n = e_n;
146 obj.du_e = kr(d1_r_u,I_v); 172
147 obj.dv_e = (obj.e_e'*Dv)'; 173 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
148 obj.du_s = (obj.e_s'*Du)'; 174 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
149 obj.dv_s = kr(I_u,d1_l_v); 175
150 obj.du_n = (obj.e_n'*Du)'; 176
151 obj.dv_n = kr(I_u,d1_r_v); 177 % Misc.
152
153 obj.x_u = x_u;
154 obj.x_v = x_v;
155 obj.y_u = y_u;
156 obj.y_v = y_v;
157
158 obj.m = m; 178 obj.m = m;
159 obj.h = [h_u h_v]; 179 obj.h = [h_u h_v];
160 obj.order = order; 180 obj.order = order;
161 obj.grid = g; 181 obj.grid = g;
162 182
163 obj.a = a; 183 obj.a = a;
164 obj.b = b; 184 obj.b = b;
165 obj.J = spdiag(J);
166 obj.Ji = spdiag(1./J);
167 obj.a11 = a11; 185 obj.a11 = a11;
168 obj.a12 = a12; 186 obj.a12 = a12;
169 obj.a22 = a22; 187 obj.a22 = a22;
170 obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv);
171 obj.lambda = lambda; 188 obj.lambda = lambda;
172
173 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
174 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
175 189
176 obj.gamm_u = h_u*ops_u.borrowing.M.d1; 190 obj.gamm_u = h_u*ops_u.borrowing.M.d1;
177 obj.gamm_v = h_v*ops_v.borrowing.M.d1; 191 obj.gamm_v = h_v*ops_v.borrowing.M.d1;
178 end 192 end
179 193