Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2d.m @ 301:d9860ebc3148 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
HypsystCurve2D Seems to work (Converges with MMS)
author Ylva Rydin <ylva.rydin@telia.com>
date Wed, 05 Oct 2016 17:36:34 +0200
parents cd30b22cee56
children 9b3d7fc61a36
comparison
equal deleted inserted replaced
300:32f3ee81bc37 301:d9860ebc3148
126 126
127 127
128 function [closure, penalty]=boundary_condition_char(obj,boundary) 128 function [closure, penalty]=boundary_condition_char(obj,boundary)
129 params=obj.params; 129 params=obj.params;
130 x=obj.x; y=obj.y; 130 x=obj.x; y=obj.y;
131 side=max(length(x),length(y)); 131
132
133 switch boundary 132 switch boundary
134 case {'w','W','west'} 133 case {'w','W','west'}
135 e_=obj.e_w; 134 e_=obj.e_w;
136 mat=obj.A; 135 mat=obj.A;
137 boundPos='l'; 136 boundPos='l';
138 Hi=obj.Hxi; 137 Hi=obj.Hxi;
139 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y); 138 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y);
139 side=max(length(y));
140 case {'e','E','east'} 140 case {'e','E','east'}
141 e_=obj.e_e; 141 e_=obj.e_e;
142 mat=obj.A; 142 mat=obj.A;
143 boundPos='r'; 143 boundPos='r';
144 Hi=obj.Hxi; 144 Hi=obj.Hxi;
145 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y); 145 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y);
146 side=max(length(y));
146 case {'s','S','south'} 147 case {'s','S','south'}
147 e_=obj.e_s; 148 e_=obj.e_s;
148 mat=obj.B; 149 mat=obj.B;
149 boundPos='l'; 150 boundPos='l';
150 Hi=obj.Hyi; 151 Hi=obj.Hyi;
151 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1)); 152 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1));
153 side=max(length(x));
152 case {'n','N','north'} 154 case {'n','N','north'}
153 e_=obj.e_n; 155 e_=obj.e_n;
154 mat=obj.B; 156 mat=obj.B;
155 boundPos='r'; 157 boundPos='r';
156 Hi=obj.Hyi; 158 Hi=obj.Hyi;
157 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end)); 159 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end));
160 side=max(length(x));
158 end 161 end
159 162
160 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 163 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
161 164
162 switch boundPos 165 switch boundPos
163 case {'l'} 166 case {'l'}
164 tau=sparse(obj.n*side,pos*side); 167 tau=sparse(obj.n*side,pos);
165 Vi_plus=Vi(1:pos*side,:); 168 Vi_plus=Vi(1:pos,:);
166 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side)); 169 tau(1:pos,:)=-abs(D(1:pos,1:pos));
167 closure=Hi*e_*V*tau*Vi_plus*e_'; 170 closure=Hi*e_*V*tau*Vi_plus*e_';
168 penalty=-Hi*e_*V*tau*Vi_plus; 171 penalty=-Hi*e_*V*tau*Vi_plus;
169 case {'r'} 172 case {'r'}
170 tau=sparse(obj.n*side,neg*side); 173 tau=sparse(obj.n*side,neg);
171 tau((pos+zeroval)*side+1:obj.n*side,:)=-abs(D((pos+zeroval)*side+1:obj.n*side,(pos+zeroval)*side+1:obj.n*side)); 174 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
172 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:); 175 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:);
173 closure=Hi*e_*V*tau*Vi_minus*e_'; 176 closure=Hi*e_*V*tau*Vi_minus*e_';
174 penalty=-Hi*e_*V*tau*Vi_minus; 177 penalty=-Hi*e_*V*tau*Vi_minus;
175 end 178 end
176 end 179 end
177 180
178 181
179 function [closure,penalty]=boundary_condition_general(obj,boundary,L) 182 function [closure,penalty]=boundary_condition_general(obj,boundary,L)
180 params=obj.params; 183 params=obj.params;
181 x=obj.x; y=obj.y; 184 x=obj.x; y=obj.y;
182 side=max(length(x),length(y));
183 185
184 switch boundary 186 switch boundary
185 case {'w','W','west'} 187 case {'w','W','west'}
186 e_=obj.e_w; 188 e_=obj.e_w;
187 mat=obj.A; 189 mat=obj.A;
188 boundPos='l'; 190 boundPos='l';
189 Hi=obj.Hxi; 191 Hi=obj.Hxi;
190 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y); 192 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y);
191 L=obj.evaluateCoefficientMatrix(L,x(1),y); 193 L=obj.evaluateCoefficientMatrix(L,x(1),y);
194 side=max(length(y));
192 case {'e','E','east'} 195 case {'e','E','east'}
193 e_=obj.e_e; 196 e_=obj.e_e;
194 mat=obj.A; 197 mat=obj.A;
195 boundPos='r'; 198 boundPos='r';
196 Hi=obj.Hxi; 199 Hi=obj.Hxi;
197 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y); 200 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y);
198 L=obj.evaluateCoefficientMatrix(L,x(end),y); 201 L=obj.evaluateCoefficientMatrix(L,x(end),y);
202 side=max(length(y));
199 case {'s','S','south'} 203 case {'s','S','south'}
200 e_=obj.e_s; 204 e_=obj.e_s;
201 mat=obj.B; 205 mat=obj.B;
202 boundPos='l'; 206 boundPos='l';
203 Hi=obj.Hyi; 207 Hi=obj.Hyi;
204 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1)); 208 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1));
205 L=obj.evaluateCoefficientMatrix(L,x,y(1)); 209 L=obj.evaluateCoefficientMatrix(L,x,y(1));
210 side=max(length(x));
206 case {'n','N','north'} 211 case {'n','N','north'}
207 e_=obj.e_n; 212 e_=obj.e_n;
208 mat=obj.B; 213 mat=obj.B;
209 boundPos='r'; 214 boundPos='r';
210 Hi=obj.Hyi; 215 Hi=obj.Hyi;
211 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end)); 216 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end));
212 L=obj.evaluateCoefficientMatrix(L,x,y(end)); 217 L=obj.evaluateCoefficientMatrix(L,x,y(end));
218 side=max(length(x));
213 end 219 end
214 220
215 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 221 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
216 222
217 switch boundPos 223 switch boundPos
218 case {'l'} 224 case {'l'}
219 tau=sparse(obj.n*side,pos*side); 225 tau=sparse(obj.n*side,pos);
220 Vi_plus=Vi(1:pos*side,:); 226 Vi_plus=Vi(1:pos,:);
221 Vi_minus=Vi(pos*side+1:obj.n*side,:); 227 Vi_minus=Vi(pos+zeroval+1:obj.n*side,:);
222 V_plus=V(:,1:pos*side); 228 V_plus=V(:,1:pos);
223 V_minus=V(:,(pos+zeroval)*side+1:obj.n*side); 229 V_minus=V(:,(pos+zeroval)+1:obj.n*side);
224 230
225 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side)); 231 tau(1:pos,:)=-abs(D(1:pos,1:pos));
226 R=-inv(L*V_plus)*(L*V_minus); 232 R=-inv(L*V_plus)*(L*V_minus);
227 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; 233 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
228 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L; 234 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L;
229 case {'r'} 235 case {'r'}
230 tau=sparse(obj.n*side,neg*side); 236 tau=sparse(obj.n*side,neg);
231 tau((pos+zeroval)*side+1:obj.n*side,:)=-abs(D((pos+zeroval)*side+1:obj.n*side,(pos+zeroval)*side+1:obj.n*side)); 237 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
232 Vi_plus=Vi(1:pos*side,:); 238 Vi_plus=Vi(1:pos,:);
233 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:); 239 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:);
234 240
235 V_plus=V(:,1:pos*side); 241 V_plus=V(:,1:pos);
236 V_minus=V(:,(pos+zeroval)*side+1:obj.n*side); 242 V_minus=V(:,(pos+zeroval)+1:obj.n*side);
237 R=-inv(L*V_minus)*(L*V_plus); 243 R=-inv(L*V_minus)*(L*V_plus);
238 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 244 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
239 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L; 245 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L;
240 end 246 end
241 end 247 end
242 248
243 249
244 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y) 250 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y)
245 params=obj.params; 251 params=obj.params;
246 syms xs ys; 252 syms xs ys
247 [V, D]=eig(mat(params,xs,ys)); 253 [V, D]=eig(mat(params,xs,ys));
248 xs=1;ys=1; 254 xs=x;
249 DD=eval(diag(D)); 255 ys=y;
250 256
251 poseig=find(DD>0);
252 zeroeig=find(DD==0);
253 negeig=find(DD<0);
254 syms xs ys
255 DD=diag(D);
256
257 D=diag([DD(poseig);DD(zeroeig); DD(negeig)]);
258 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)];
259 xs=x; ys=y;
260
261 side=max(length(x),length(y)); 257 side=max(length(x),length(y));
262 Dret=zeros(obj.n,side*obj.n); 258 Dret=zeros(obj.n,side*obj.n);
263 Vret=zeros(obj.n,side*obj.n); 259 Vret=zeros(obj.n,side*obj.n);
264 for ii=1:obj.n 260 for ii=1:obj.n
265 for jj=1:obj.n 261 for jj=1:obj.n
267 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii)); 263 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii));
268 end 264 end
269 end 265 end
270 266
271 D=sparse(Dret); 267 D=sparse(Dret);
272 V=sparse(normc(Vret)); 268 V=sparse(Vret);
273 V=obj.evaluateCoefficientMatrix(V,x,y); 269 V=obj.evaluateCoefficientMatrix(V,x,y);
274 D=obj.evaluateCoefficientMatrix(D,x,y); 270 D=obj.evaluateCoefficientMatrix(D,x,y);
271 DD=diag(D);
272
273 poseig=(DD>0);
274 zeroeig=(DD==0);
275 negeig=(DD<0);
276
277 D=diag([DD(poseig); DD(zeroeig); DD(negeig)]);
278 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)];
275 Vi=inv(V); 279 Vi=inv(V);
276 signVec=[length(poseig),length(zeroeig),length(negeig)]; 280 signVec=[sum(poseig),sum(zeroeig),sum(negeig)];
277 end 281 end
278 282
279 end 283 end
280 end 284 end