Mercurial > repos > public > sbplib

annotate +scheme/Hypsyst2d.m @ 1333:0aefcb30cab4 feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add support for RK6
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sat, 07 May 2022 10:30:59 +0200
parents 8d73fcdb07a5
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
1 classdef Hypsyst2d < scheme.Scheme
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
2 properties
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
3 m % Number of points in each direction, possibly a vector
293
2d604d16842c Works with varying coefficients and char boundary condition
Ylva Rydin <ylva.rydin@telia.com>
parents: 292
diff changeset
4 n %size of system
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
5 h % Grid spacing
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
6 x,y % Grid
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
7 X,Y % Values of x and y for each grid point
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
8 order % Order accuracy for the approximation
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
9
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
10 D % non-stabalized scheme operator
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
11 A, B, E %Coefficient matrices
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
12
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
13 H % Discrete norm
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
14 % Norms in the x and y directions
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
15 Hxi,Hyi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
16 I_x,I_y, I_N
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
17 e_w, e_e, e_s, e_n
297
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
18 params %parameters for the coeficient matrice
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
19 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
20
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
21 methods
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
22 %Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
23 function obj = Hypsyst2d(m, lim, order, A, B, E, params)
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
24 default_arg('E', [])
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
25 xlim = lim{1};
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
26 ylim = lim{2};
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
27
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
28 if length(m) == 1
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
29 m = [m m];
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
30 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
31
297
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
32 obj.A=A;
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
33 obj.B=B;
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
34 obj.E=E;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
35
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
36 m_x = m(1);
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
37 m_y = m(2);
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
38 obj.params = params;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
39
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
40 ops_x = sbp.D2Standard(m_x,xlim,order);
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
41 ops_y = sbp.D2Standard(m_y,ylim,order);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
42
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
43 obj.x = ops_x.x;
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
44 obj.y = ops_y.x;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
45
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
46 obj.X = kr(obj.x,ones(m_y,1));
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
47 obj.Y = kr(ones(m_x,1),obj.y);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
48
297
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
49 Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y);
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
50 Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y);
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
51 Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
52
297
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
53 obj.n = length(A(obj.params,0,0));
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
54
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
55 I_n = eye(obj.n);I_x = speye(m_x);
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
56 obj.I_x = I_x;
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
57 I_y = speye(m_y);
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
58 obj.I_y = I_y;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
59
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
60
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
61 D1_x = kr(I_n, ops_x.D1, I_y);
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
62 obj.Hxi = kr(I_n, ops_x.HI, I_y);
297
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
63 D1_y = kr(I_n, I_x, ops_y.D1);
cd30b22cee56 Have tried to make a curvelinear sheme for hypsysts. Does not really work yet...
Ylva Rydin <ylva.rydin@telia.com>
parents: 296
diff changeset
64 obj.Hyi = kr(I_n, I_x, ops_y.HI);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
65
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
66 obj.e_w = kr(I_n, ops_x.e_l, I_y);
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
67 obj.e_e = kr(I_n, ops_x.e_r, I_y);
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
68 obj.e_s = kr(I_n, I_x, ops_y.e_l);
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
69 obj.e_n = kr(I_n, I_x, ops_y.e_r);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
70
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
71 obj.m = m;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
72 obj.h = [ops_x.h ops_y.h];
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
73 obj.order = order;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
74
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
75 obj.D = -Aevaluated*D1_x-Bevaluated*D1_y-Eevaluated;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
76
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
77 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
78
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
79 % Closure functions return the opertors applied to the own doamin to close the boundary
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
80 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
81 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
82 % type is a string specifying the type of boundary condition if there are several.
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
83 % data is a function returning the data that should be applied at the boundary.
293
2d604d16842c Works with varying coefficients and char boundary condition
Ylva Rydin <ylva.rydin@telia.com>
parents: 292
diff changeset
84 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
85 default_arg('type','char');
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
86 switch type
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
87 case{'c','char'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
88 [closure,penalty] = boundary_condition_char(obj,boundary);
293
2d604d16842c Works with varying coefficients and char boundary condition
Ylva Rydin <ylva.rydin@telia.com>
parents: 292
diff changeset
89 case{'general'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
90 [closure,penalty] = boundary_condition_general(obj,boundary,L);
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
91 otherwise
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
92 error('No such boundary condition')
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
93 end
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
94 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
95
946
706d1c2b4199 Raname opts to type in a bunch of interface methods
Jonatan Werpers <jonatan@werpers.com>
parents: 910
diff changeset
96 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
706d1c2b4199 Raname opts to type in a bunch of interface methods
Jonatan Werpers <jonatan@werpers.com>
parents: 910
diff changeset
97 error('Not implemented');
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
98 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
99
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
100 function N = size(obj)
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
101 N = obj.m;
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
102 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
103
295
da0131655035 Fixed some formatting and naming.
Jonatan Werpers <jonatan@werpers.com>
parents: 294
diff changeset
104 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y)
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
105 params = obj.params;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
106
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
107 if isa(mat,'function_handle')
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
108 [rows,cols] = size(mat(params,0,0));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
109 matVec = mat(params,X',Y');
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
110 matVec = sparse(matVec);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
111 side = max(length(X),length(Y));
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
112 else
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
113 matVec = mat;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
114 [rows,cols] = size(matVec);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
115 side = max(length(X),length(Y));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
116 cols = cols/side;
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
117 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
118 ret = cell(rows,cols);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
119
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
120 for ii = 1:rows
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
121 for jj=1:cols
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
122 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side));
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
123 end
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
124 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
125 ret = cell2mat(ret);
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
126 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
127
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
128 %Characteristic boundary conditions
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
129 function [closure, penalty] = boundary_condition_char(obj,boundary)
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
130 params = obj.params;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
131 x = obj.x;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
132 y = obj.y;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
133
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
134 switch boundary
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
135 case {'w','W','west'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
136 e_ = obj.e_w;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
137 mat = obj.A;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
138 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
139 Hi = obj.Hxi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
140 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(1),y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
141 side = max(length(y));
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
142 case {'e','E','east'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
143 e_ = obj.e_e;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
144 mat = obj.A;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
145 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
146 Hi = obj.Hxi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
147 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(end),y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
148 side = max(length(y));
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
149 case {'s','S','south'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
150 e_ = obj.e_s;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
151 mat = obj.B;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
152 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
153 Hi = obj.Hyi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
154 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(1));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
155 side = max(length(x));
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
156 case {'n','N','north'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
157 e_ = obj.e_n;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
158 mat = obj.B;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
159 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
160 Hi = obj.Hyi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
161 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
162 side = max(length(x));
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
163 end
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
164 pos = signVec(1);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
165 zeroval = signVec(2);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
166 neg = signVec(3);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
167
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
168 switch boundPos
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
169 case {'l'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
170 tau = sparse(obj.n*side,pos);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
171 Vi_plus = Vi(1:pos,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
172 tau(1:pos,:) = -abs(D(1:pos,1:pos));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
173 closure = Hi*e_*V*tau*Vi_plus*e_';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
174 penalty = -Hi*e_*V*tau*Vi_plus;
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
175 case {'r'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
176 tau = sparse(obj.n*side,neg);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
177 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
178 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
179 closure = Hi*e_*V*tau*Vi_minus*e_';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
180 penalty = -Hi*e_*V*tau*Vi_minus;
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
181 end
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
182 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
183
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
184 % General boundary condition in the form Lu=g(x)
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
185 function [closure,penalty] = boundary_condition_general(obj,boundary,L)
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
186 params = obj.params;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
187 x = obj.x;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
188 y = obj.y;
997
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
189 e_ = obj.getBoundaryOperator('e', boundary);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
190
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
191 switch boundary
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
192 case {'w','W','west'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
193 mat = obj.A;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
194 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
195 Hi = obj.Hxi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
196 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(1),y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
197 L = obj.evaluateCoefficientMatrix(L,x(1),y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
198 side = max(length(y));
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
199 case {'e','E','east'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
200 mat = obj.A;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
201 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
202 Hi = obj.Hxi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
203 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(end),y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
204 L = obj.evaluateCoefficientMatrix(L,x(end),y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
205 side = max(length(y));
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
206 case {'s','S','south'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
207 mat = obj.B;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
208 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
209 Hi = obj.Hyi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
210 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(1));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
211 L = obj.evaluateCoefficientMatrix(L,x,y(1));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
212 side = max(length(x));
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
213 case {'n','N','north'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
214 mat = obj.B;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
215 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
216 Hi = obj.Hyi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
217 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end));
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
218 L = obj.evaluateCoefficientMatrix(L,x,y(end));
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
219 side = max(length(x));
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
220 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
221
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
222 pos = signVec(1);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
223 zeroval = signVec(2);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
224 neg = signVec(3);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
225
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
226 switch boundPos
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
227 case {'l'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
228 tau = sparse(obj.n*side,pos);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
229 Vi_plus = Vi(1:pos,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
230 Vi_minus = Vi(pos+zeroval+1:obj.n*side,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
231 V_plus = V(:,1:pos);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
232 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
233
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
234 tau(1:pos,:) = -abs(D(1:pos,1:pos));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
235 R = -inv(L*V_plus)*(L*V_minus);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
236 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
237 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L;
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
238 case {'r'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
239 tau = sparse(obj.n*side,neg);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
240 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
241 Vi_plus = Vi(1:pos,:);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
242 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
243
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
244 V_plus = V(:,1:pos);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
245 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
246 R = -inv(L*V_minus)*(L*V_plus);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
247 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
248 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
294
8ff6ec6249e8 "General" boundary conditions implemented
Ylva Rydin <ylva.rydin@telia.com>
parents: 293
diff changeset
249 end
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
250 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
251
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
252 % Function that diagonalizes a symbolic matrix A as A=V*D*Vi
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
253 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
254 % [d+ ]
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
255 % D = [ d0 ]
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
256 % [ d-]
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
257 % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
258 function [V,Vi, D,signVec] = matrixDiag(obj,mat,x,y)
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
259 params = obj.params;
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 297
diff changeset
260 syms xs ys
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
261 [V, D]= eig(mat(params,xs,ys));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
262 Vi = inv(V);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
263 xs = x;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
264 ys = y;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
265
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
266 side = max(length(x),length(y));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
267 Dret = zeros(obj.n,side*obj.n);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
268 Vret = zeros(obj.n,side*obj.n);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
269 Viret = zeros(obj.n,side*obj.n);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
270
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
271 for ii = 1:obj.n
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
272 for jj = 1:obj.n
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
273 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
274 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii));
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
275 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
276 end
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
277 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
278
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
279 D = sparse(Dret);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
280 V = sparse(Vret);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
281 Vi = sparse(Viret);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
282 V = obj.evaluateCoefficientMatrix(V,x,y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
283 Vi = obj.evaluateCoefficientMatrix(Vi,x,y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
284 D = obj.evaluateCoefficientMatrix(D,x,y);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
285 DD = diag(D);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
286
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
287 poseig = (DD>0);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
288 zeroeig = (DD==0);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
289 negeig = (DD<0);
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
290
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
291 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
292 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
293 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
294 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
291
807dfe8be3ec Have made a lot of stupid changes in hypsyst in order to find a stupid bug
Ylva Rydin <ylva.rydin@telia.com>
parents: 290
diff changeset
295 end
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
296
997
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
297 % Returns the boundary operator op for the boundary specified by the string boundary.
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
298 % op -- string or a cell array of strings
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
299 % boundary -- string
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
300 function varargout = getBoundaryOperator(obj, op, boundary)
1042
8d73fcdb07a5 Add asserts to boundary identifier inputs
Jonatan Werpers <jonatan@werpers.com>
parents: 997
diff changeset
301 assertIsMember(boundary, {'w', 'e', 's', 'n'})
997
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
302
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
303 if ~iscell(op)
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
304 op = {op};
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
305 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
306
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
307 for i = 1:numel(op)
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
308 switch op{i}
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
309 case 'e'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
310 switch boundary
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
311 case 'w'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
312 e = obj.e_w;
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
313 case 'e'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
314 e = obj.e_e;
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
315 case 's'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
316 e = obj.e_s;
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
317 case 'n'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
318 e = obj.e_n;
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
319 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
320 varargout{i} = e;
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
321 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
322 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
323 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
324
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
325 % Returns square boundary quadrature matrix, of dimension
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
326 % corresponding to the number of boundary points
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
327 %
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
328 % boundary -- string
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
329 function H_b = getBoundaryQuadrature(obj, boundary)
1042
8d73fcdb07a5 Add asserts to boundary identifier inputs
Jonatan Werpers <jonatan@werpers.com>
parents: 997
diff changeset
330 assertIsMember(boundary, {'w', 'e', 's', 'n'})
997
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
331
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
332 e = obj.getBoundaryOperator('e', boundary);
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
333
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
334 switch boundary
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
335 case 'w'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
336 H_b = inv(e'*obj.Hyi*e);
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
337 case 'e'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
338 H_b = inv(e'*obj.Hyi*e);
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
339 case 's'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
340 H_b = inv(e'*obj.Hxi*e);
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
341 case 'n'
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
342 H_b = inv(e'*obj.Hxi*e);
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
343 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
344 end
78db023a7fe3 Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
Martin Almquist <malmquist@stanford.edu>
parents: 946
diff changeset
345
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
346 end
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
347 end