Mercurial > repos > public > sbplib

annotate +scheme/Hypsyst2d.m @ 1302:a0d615bde7f8 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add the hollow option to the anisotropic diffops
author Martin Almquist <malmquist@stanford.edu>
date Fri, 10 Jul 2020 20:24:23 -0700
parents 706d1c2b4199
children 78db023a7fe3
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;
905
459eeb99130f Include type as (optional) input parameter in the interface method of all schemes.
Martin Almquist <malmquist@stanford.edu>
parents: 369
diff changeset
189
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
190 switch boundary
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
191 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
192 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
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 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
201 mat = obj.A;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
202 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
203 Hi = obj.Hxi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
204 [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
205 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
206 side = max(length(y));
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
207 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
208 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
209 mat = obj.B;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
210 boundPos = 'l';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
211 Hi = obj.Hyi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
212 [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
213 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
214 side = max(length(x));
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
215 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
216 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
217 mat = obj.B;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
218 boundPos = 'r';
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
219 Hi = obj.Hyi;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
220 [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
221 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
222 side = max(length(x));
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
223 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
224
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
225 pos = signVec(1);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
226 zeroval = signVec(2);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
227 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
228
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
229 switch boundPos
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
230 case {'l'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
231 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
232 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
233 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
234 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
235 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
236
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
237 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
238 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
239 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
240 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
241 case {'r'}
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
242 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
243 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
244 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
245 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
246
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
247 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
248 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
249 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
250 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
251 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
252 end
292
3d275c5e45b3 Changed how the matrices are built
Ylva Rydin <ylva.rydin@telia.com>
parents: 291
diff changeset
253 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
254
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
255 % 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
256 % 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
257 % [d+ ]
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
258 % D = [ d0 ]
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
259 % [ 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
260 % 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
261 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
262 params = obj.params;
301
d9860ebc3148 HypsystCurve2D Seems to work (Converges with MMS)
Ylva Rydin <ylva.rydin@telia.com>
parents: 297
diff changeset
263 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
264 [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
265 Vi = inv(V);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
266 xs = x;
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
267 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
268
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
269 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
270 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
271 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
272 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
273
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
274 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
275 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
276 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
277 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
278 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
279 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
280 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
281
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
282 D = sparse(Dret);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
283 V = sparse(Vret);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
284 Vi = sparse(Viret);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
285 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
286 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
287 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
288 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
289
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
290 poseig = (DD>0);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
291 zeroeig = (DD==0);
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
292 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
293
369
9d1fc984f40d Added some comments and cleaned up the code a little
Ylva Rydin <ylva.rydin@telia.com>
parents: 352
diff changeset
294 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
295 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
296 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
297 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
298 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
299
290
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
300 end
d32f674bcbe5 A first attempt to make a general scheme fo hyperbolic systems
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
301 end